Hard and soft materials: putting consistent van der Waals density functionals to work

The vdW-DF method is a systematic approach to design XC functionals that capture truly nonlocal correlation effects. Pure vdW interactions (produced by electrodynamic coupling of electron–hole pairs [25, 31, 39, 71, 72]) are examples of nonlocal-correlation effects. Another example is the screening (by itinerant valence electrons) that shifts orbital energies as, for example, captured in a cumulant expansion [73]. In the vdW-DF design we note that both are reflected in an electrodynamics reformulation of the XC functional [25]. This allows us to treat all XC effects on the same footing in the electron-gas tradition.

In practice, we use a GGA-type functional ${E}_{\mathrm{XC}}^{\text{in}}$ to define an effective (nonlocal) model of the frequency-dependence of the electron-gas susceptibility α(ω). For reasons discussed elsewhere [22, 36, 39], we limit this input to local-density approximation (LDA) plus a simple approximation for gradient-corrected exchange. We formally express the internal semilocal functional ${E}_{\mathrm{XC}}^{\text{in}}$ as the trace of a plasmon propagator S_XC( r , r', ω),

$$\begin{equation}{E}_{\mathrm{XC}}^{\text{in}}={\int }_{0}^{\infty }\frac{\mathrm{d}u}{2\pi }\enspace \mathrm{Tr}\left\{{S}_{\mathrm{XC}}(\omega =\text{i}u)\right\}-{E}_{\text{self}}.\end{equation} \tag{ 1 }$$

The trace is here taken over the spatial coordinates of S_XC( r , r'; ω). The term E_self denotes an infinite self-energy that removes the formal divergence. For spin-carrying systems we work with the spin-up and spin-down density components n_s=↑,↓( r ) of the total electron density n( r ) = n_↑( r ) + n_↓( r ). The spin polarization η( r ) = [n_↑( r ) − n_↓( r )]/[n_↑( r ) + n_↓( r )] impacts the GGA-type internal functional ${E}_{\mathrm{XC}}^{\text{in}}$, and must therefore be directly reflected in the details of the plasmon propagator S_XC( r , r', ω) [38].

The key point is that the model plasmon propagator S_XC also defines an effective GGA-level model dielectric function (iu) = exp(S_XC(iu)) and a corresponding model susceptibility α(iu) = ((iu) − 1)/4π. Moreover, by enforcing current conservation, the dielectrics modeling also defines the full vdW-DF specification of the XC functional,

$$\begin{equation}{E}_{\mathrm{XC}}^{\text{DF}}={\int }_{0}^{\infty }\frac{\mathrm{d}u}{2\pi }\enspace \mathrm{Tr}\left\{\mathrm{ln}(\nabla {\epsilon}(\text{i}u)\cdot \nabla G)\right\}-{E}_{\text{self}},\end{equation} \tag{ 2 }$$

where G denotes the Coulomb Green function. By expanding equation (2) to first order in S_XC, one formally recoups the internal GGA-type functional equation (1). By further expanding to second order, we obtain the vdW-DF determination of corresponding nonlocal-correlation effects,

$$\begin{equation}{E}_{\text{c}}^{\text{nl,}\;\text{sp}}={\int }_{0}^{\infty }\frac{\mathrm{d}u}{4\pi }\enspace \mathrm{Tr}\left\{{S}_{\mathrm{XC}}^{2}-{(\nabla {S}_{\mathrm{XC}}\cdot \nabla G)}^{2}\right\}.\end{equation} \tag{ 3 }$$

As indicated by superscript 'sp', the nonlocal-correlation term depends on the spatial variation in the spin polarization η( r ) through S_XC( r , r'; ω) [38].

Functionals of the vdW-DF family are generally expressed as

$$\begin{equation}{E}_{\mathrm{XC}}^{\text{vdW}-\text{DF\#}}={E}_{\mathrm{XC}}^{0}+{E}_{\text{c}}^{\text{nl,}\;\text{sp}},\end{equation} \tag{ 4 }$$

where ${E}_{\mathrm{XC}}^{0}={E}_{\mathrm{XC}}^{\text{in}}+\delta {E}_{\text{x}}^{0}$ contains nothing but LDA and the gradient-corrected exchange while ${E}_{\text{c}}^{\text{nl,}\;\text{sp}}$ is the nonlocal-correlation term. The Lindhard–Dyson screening logic formally mandates that the cross-over exchange term $\delta {E}_{\text{x}}^{0}$ must vanish, thus setting the balance between exchange and correlation [25]. There are, however, practical limitations that prevent us from going directly for such fully consistent implementations. In the consistent-exchange vdW-DF-cx version we have chosen $\delta {E}_{\text{x}}^{0}$ so that the nonzero cross-over term does not affect binding energies in typical bulk and in typical molecular-interaction cases, as discussed separately in references [4, 25, 26, 74].

For actual evaluations, we use a two-pole approximation for the shape of the plasmon-pole propagator S_XC. This plasmon-pole description, and hence the resulting vdW-DF version, is effectively set by the choice of the semilocal internal functional E_XC via equation (1). This leads to the computationally efficient nonlocal-correlation determination

$$\begin{equation}{E}_{\text{c}}^{\text{nl}}=\frac{1}{2}{\int }_{\boldsymbol{r}}{\int }_{{\boldsymbol{r}}^{\prime }}n(\boldsymbol{r})n({\boldsymbol{r}}^{\prime }){\Phi}(D;{q}_{0}(\boldsymbol{r});{q}_{0}({\boldsymbol{r}}^{\prime })),\end{equation} \tag{ 5 }$$

where D = | r − r'|. It is given by a universal kernel form Φ, as discussed in reference [75]. In equation (5), the values of q₀( r ) and q₀( r') characterize the model plasmon dispersion.

The above-summarized vdW-DF framework leaves no ambiguity about how to incorporate spin-polarization effects in ${E}_{\text{c}}^{\text{nl,}\;\text{sp}}$. Spin enters via the exchange and via the LDA-correlation parts of ${E}_{\mathrm{XC}}^{\text{in}}$, as given by the spin-scaling description and by the now-standard PW92 formulation of LDA [76]. This, in turn, determines the form of S_XC and ultimately equation (5). Specifically, the values of q₀( r ) and q₀( r') are set from the local energy density of the internal semi-local functional to reflect the density and spin impact on the underlying screening description. Of course, it is imperative to also include spin in the ${E}_{\mathrm{XC}}^{0}$ term. This proper spin vdW-DF formulation [38] is implemented in QUANTUM ESPRESSO [41, 42] and in the LIBVDWXC library [70], permitting fixed-cell calculations but not (until now) general variable-cell vdW-DF studies.

For non-spin-polarized problems, there has since long existed a formulation of stress in vdW-DF calculations [77]. This allows effective KS structure optimizations as implemented in QUANTUM ESPRESSO. We now present a spin vdW-DF extension of stress to enhance the KS-structure search part. It is based on the ideas of Nielsen and Martin [78] and we first summarize the non-spin vdW-DF stress calculations, as derived by Sabatini and co-workers [77]. We consider the impact of unit-cell and coordinate scaling, for example, as expressed in Cartesian coordinates for a position vector ${r}_{\alpha }\to {\tilde{r}}_{\alpha }={\sum }_{\beta }({\delta }_{\alpha ,\beta }+{\varepsilon }_{\alpha ,\beta }){r}_{\beta }$, where ɛ_α,β is the strain tensor. This scaling affects the double Jacobian, the total-electron density n( r'), the total-density gradient ∇n( r ) and the coordinate-separation variable D inside Φ in equation (5). Details of these different scaling effects are discussed in appendix A for the spin-polarized case.

In the absence of spin polarization, Sabatini and co-workers [77] derived the nonlocal-correlation stress tensor contribution

$$\begin{align}\hfill {\sigma }_{\text{c},\alpha ,\beta }^{\text{nl}}& ={\delta }_{\alpha ,\beta }\left[2{E}_{\text{c}}^{\text{nl}}-{\int }_{\boldsymbol{r}}n(\boldsymbol{r}){v}_{\text{c}}^{\text{nl}}(\boldsymbol{r})\right]+\frac{1}{2}{\int }_{\boldsymbol{r}}{\int }_{{\boldsymbol{r}}^{\prime }}n(\boldsymbol{r})n({\boldsymbol{r}}^{\prime })\frac{\partial {\Phi}}{\partial D}{C}_{\alpha ,\beta }(\boldsymbol{r},{\boldsymbol{r}}^{\prime })\hfill \\ \hfill & \quad -{\int }_{\boldsymbol{r}}{\int }_{{\boldsymbol{r}}^{\prime }}n(\boldsymbol{r})n({\boldsymbol{r}}^{\prime })\frac{\partial {\Phi}}{\partial {q}_{0}}{G}_{\alpha ,\beta }(\boldsymbol{r}),\hfill \end{align} \tag{ 6 }$$

where

$$\begin{equation}{C}_{\alpha ,\beta }(\boldsymbol{r},{\boldsymbol{r}}^{\prime })=({r}_{\alpha }-{r}_{\alpha }^{\prime })({r}_{\beta }-{r}_{\beta }^{\prime })/D,\end{equation} \tag{ 7 }$$

and

$$\begin{equation}{G}_{\alpha ,\beta }(\boldsymbol{r})=\frac{\partial {q}_{0}(\boldsymbol{r})}{\partial \vert \nabla n(\boldsymbol{r})\vert }\frac{(\partial n(\boldsymbol{r})/\partial {r}_{\alpha })(\partial n(\boldsymbol{r})/\partial {r}_{\beta })}{\vert \nabla n(\boldsymbol{r})\vert }.\end{equation} \tag{ 8 }$$

This stress component is in part given by the nonlocal-correlation contributions, ${v}_{\text{c}}^{\text{nl}}(\boldsymbol{r})$ and ${E}_{\text{c}}^{\text{nl}}$, to the XC potential and XC energy. These contributions are given by local values of an inverse length scale, q₀( r ), that determines the local plasmon dispersion. As such, the contributions depend on the density gradients and hence have an indirect dependence on coordinate scaling, as summarized in equation (8).

For stress evaluations it is important to note that the density gradient will scale both since the density scales with the unit-cell size and because the formal expression for the spatial gradient scales with coordinate dilation even at a fixed density. The former effect is incorporated by the term containing ${v}_{\text{c}}^{\text{nl}}(\boldsymbol{r})$. The latter effect is captured by the term in the last row of equation (6). Fortunately, the local variation of q₀ is already computed as part of any self-consistent (spin) vdW-DF calculation.

The internal functional depends on the variation in the spin polarization η( r ) and this spin dependence impacts the plasmon dispersion, and ultimately ${E}_{\text{c}}^{\text{nl,}\;\text{sp}}$, through the local q₀( r ) values. To also compute stresses in spin vdW-DF we must update equation (6) accordingly; details are discussed in appendix A.

We find that the second line of equation (6) only changes to the extent that the values of the q₀'s must now be evaluated for η( r ) ≠ 0. For the first line, there is a small modification since separate XC potentials now act on n_↑( r ) and n_↓( r ).

Finally, for an update of the third line of equation (6), we simply track the variation of the local q₀ values on both spin-density gradient terms. Thus, the resulting spin-vdW-DF stress tensor expression becomes

$$\begin{align}\hfill {\sigma }_{\text{c},\alpha ,\beta }^{\text{nl,}\;\text{sp}}& ={\delta }_{\alpha ,\beta }\left[2{E}_{\text{c}}^{\text{nl}}-\sum\limits _{\text{s}={\uparrow},{\downarrow}}{\int }_{\boldsymbol{r}}{n}_{\text{s}}(\boldsymbol{r}){v}_{\text{c,s}}^{\text{nl,}\;\text{sp}}(\boldsymbol{r})\right]+\frac{1}{2}{\int }_{\boldsymbol{r}}{\int }_{{\boldsymbol{r}}^{\prime }}n(\boldsymbol{r})n({\boldsymbol{r}}^{\prime })\frac{\partial {\Phi}}{\partial D}{C}_{\alpha ,\beta }(\boldsymbol{r},{\boldsymbol{r}}^{\prime })\hfill \\ \hfill & \quad -{\int }_{\boldsymbol{r}}{\int }_{{\boldsymbol{r}}^{\prime }}n(\boldsymbol{r})n({\boldsymbol{r}}^{\prime })\frac{\partial {\Phi}}{\partial {q}_{0}}\sum\limits _{\text{s}={\uparrow},{\downarrow}}{G}_{\alpha ,\beta }^{\text{s}}(\boldsymbol{r}),\hfill \end{align} \tag{ 9 }$$

where

$$\begin{equation}{G}_{\alpha ,\beta }^{\text{s}={\uparrow},{\downarrow}}(\boldsymbol{r})=\frac{\partial {q}_{0}(\boldsymbol{r})}{\partial \vert \nabla {n}_{\text{s}}(\boldsymbol{r})\vert }\frac{(\partial {n}_{\text{s}}/\partial {r}_{\alpha })(\partial {n}_{\text{s}}/\partial {r}_{\beta })}{\vert \nabla n\vert }.\end{equation} \tag{ 10 }$$

As in the non-spin-polarized case, the spin vdW-DF stress contribution, equation (10), is conveniently given by quantities that are already computed in a self-consistent determination of the electron density variation in spin vdW-DF.

Sabatini et al [77] incorporated their non-spin stress result in the QUANTUM ESPRESSO DFT code suite. We have now completed a full spin vdW-DF implementation by having already contributed our spin-vdW-DF-stress code extension to this open-source DFT suite. Our public release makes a previous work-around (available through a compiler flag) obsolete in QUANTUM ESPRESSO: users no longer have to omit the spin impact on the evaluation of the nonlocal-correlation energy term when performing variable-cell vdW-DF calculations in spin-polarized systems.

We use and compare results obtained in PBE- and CX-associated hybrids, both unscreened and in an RSH form. We see these functionals as tool chains of closely related functionals, having either a semilocal or a truly-nonlocal correlation term. The former tool chain comprises PBE [33], PBE0 [9, 11], and HSE [12]. The latter tool chain is defined by consistent vdW-DFs [25] and is new. It comprises CX [4, 26], CX0P [14], and AHCX [15].

The PBE0 [9, 11, 13] is an unscreened hybrid based on PBE. One merely replaces 25% of the PBE exchange component with an evaluation based on the KS orbitals. The HSE functional [12] is an RSH extension of the PBE functional. We use it with 25% Fock exchange and a range separation that is described by a screening parameter μ = 0.2 Å⁻¹. This parameter defines an error-function weighting erf(μr)/r of the Coulomb interaction [12], thus limiting the Fock-exchange inclusion to short separations.

The vdW-DF-cx0 is an unscreened hybrid class that is formulated in analogy with PBE0 [9, 11, 13] but starting instead with CX. We use the zero-parameter form CX0P [14] in which the extent of Fock exchange mixing is kept fixed at 20%, following an analysis of the CX coupling-constant scaling [40] that enters the general hybrid design logic [11, 14].

Finally, the AHCX is the recently launched CX-based RSH [15]. It resembles the CX0P in that we keep the Fock exchange fraction fixed at 20% and it resembles HSE in that we keep the screening parameter fixed at the standard HSE value [12]. The screening makes AHCX calculations relevant for metallic systems [15].

For the cubic forms, figure 3 (stable above 105 K in SrTiO₃ and in general for BaZrO₃) we determine the structure and characterize the quadratic variation of the energy with atomic deformations for a range of functionals. We focus on the PBE/HSE part of the semilocal-correlation tool chain and the CX/CX0P part of the nonlocal-correlation tool chain, but also include LDA for a reference. The computed data allow us to in turn make functional-specific predictions of SrTiO₃ and BaZrO₃ properties, such as the dielectric constant (following reference [55]) and both the Γ₁₅ and R₂₅ modes. The former reflects a possibility for a ferroelectric transition while the latter reflects a potential AFD transition. Our survey and comparisons furthermore include calculations of the BaZrO₃ and SrTiO₃ elastic coefficients.

For our calculations, we use the PAW method and the VASP [59, 60] software. Our BaZrO₃ and SrTiO₃ studies are converged with respect to the wavefunction energy cutoff. We have previously shown [55] that the AFD mode is very sensitive to the oxygen PAW potential and the energy cutoff, and thus we use the hard setup in VASP. For BaZrO₃ we use an energy cutoff of 1200 eV for LDA and GGA and 1600 eV for HSE, CX, and CX0P. For SrTiO₃ energy cutoffs at 1600 eV are used for all functionals. Convergence turns out to be less sensitive to the k-point sampling and a 6 × 6 × 6 Monkhorst–Pack [103] grid was deemed sufficient for the hybrid functionals, while 8 × 8 × 8 was used for non-hybrid studies.

We also use finite-difference phonopy [116] to determine the vibrational modes and frequencies in the harmonic approximation that forms the starting point for our discussion of the cubic-SrTiO₃ and cubic-BaZrO₃ properties. We compute the R₂₅ and Γ₁₅ frequencies using the frozen phonon method with the default displacement of 0.01 Å in a 40 atom cell (being a 2 × 2 × 2 repetition of the basic cell).

The present results extend a previous BaZrO₃-only characterization [54, 55], by seeking a simple modeling-based approach to connect DFT calculations (obtained at the BO lattice constants) with measurements. The modeling is important since the experiments to which we compare are generally obtained at room or at least elevated temperatures for SrTiO₃ (which does not remain cubic below 105 K); the measurements will, in any case, be affected by the expansion that is caused by zero-point vibrational effects [55]. The modeling is, however, difficult because of the SrTiO₃ phase transition: we cannot simply track the thermal impact on structure like in references [51, 54, 55] and we cannot here make a direct comparison of thermal-expansion results among functionals or between the two cubic perovskites. Instead, we rely on extrapolation in combination with an assessment of how the vibrational ZPE impacts a cubic-perovskite structure.

An experiment-based estimate for the zero-temperature value of the lattice constant of (metastable) cubic SrTiO₃ is a_T→0 = 3.894 Å [108] (as listed in table 3). The estimate is based on tracking the temperature variation from the room temperature value a₂₉₈ = 3.905 Å [108, 123] and down to the value 3.898 Å at the SrTiO₃ phase transition temperature 105 K [107, 108]. Meanwhile, we can use our previous BaZrO₃ experience [54, 55] to extract an estimate of the vibrational ZPE effects: from the BaZrO₃ calculation we expect that the SrTiO₃ lattice constants will expand (off of a given BO result) by about 0.008 Å at T → 0 [55]. In other word, for SrTiO₃ we arrive at a back-corrected experimental lattice-constant value ${a}_{\text{ref}}^{\ast }=3.886$ Å. We can use this value to directly benchmark the accuracy of the set of functional-specific BO (cubic-)structure determinations, a₀ (listed in table 3).

Table 3. Comparison of equilibrium lattice constant a₀ (as obtained from a Birch–Murnaghan fit to DFT calculations in the BO approximation), frequencies of the lowest (and possibly soft) mode at the Γ- and R-point, Landau-model expansion parameters ($\bar{\kappa },\bar{\alpha }$, see equation (B3)), characteristic life-time (τ_R ) of the phonon trapped in the double well potential, elastic coefficients, and the high-frequency (component of the) dielectric constant ɛ_∞. We list results computed in regular and hybrid functionals having either a semilocal- (PBE and HSE) or a truly nonlocal-correlation (CX and CX0P) description; LDA results and experimental values are included for reference. All results are obtained for cubic cells at the functional-specific BO lattice constant a₀, except $\bar{\kappa }$, $\bar{\alpha }$ and τ_R which are obtained from the set of PES results, i.e., computed by deforming a 40-atom unit cell. Imaginary values in Γ₁₅ or R₂₅ reflect a possible or incipient instability in the description of the cubic cell by that functional. Corresponding identifier 'soft' in the experimental column reflects the observation of a low-temperature phase transition.

	Exper.	LDA	PBE	CX	HSE	CX0P
BaZrO3
a0 (Å)	4.188 a	4.160	4.237	4.200	4.200	4.183
Γ15 (meV)	15.2 b	13.04	11.97	13.70	13.47	14.85
R25 (meV)	5.9 a	i7.17	2.30	i2.53	6.12	5.49
$\bar{\kappa }$ (meV Å− 2 u−1)	8.51 a	−12.5	1.27	−1.30	8.54	6.68
$\bar{\alpha }$ (meV Å−4 u−2)	—	2.55	2.03	2.19	2.32	2.3
τR (10−12 s)	—	6.0	—	—	—	—
C11 (GPa)	282 c /332 d	348	290	324	298	340
C12 (GPa)	88 c	88	79	89	84	87
C44 (GPa)	97 c /97 d	90	85	88	94	97
ɛ∞	4.928 e	4.92	4.88	4.88	4.25	4.32
SrTiO3
	Exper.	LDA	PBE	CX	HSE	CX0P
a0 (Å)	3.894 f	3.862	3.942	3.905	3.900	3.880
Γ15 (meV)	soft	7.16	i17.34	i6.69	i15.52	i1.38
R25 (meV)	soft	i11.15	i8.34	i8.88	i3.29	i5.56
$\bar{\kappa }$ (meV Å− 2 u−1)	soft	−29.2	−16.3	−18.4	−3.8	−8.9
$\bar{\alpha }$ (meV Å− 4 u−1)	—	3.36	2.84	3.07	3.17	3.34
τR (10−12 s)	—	14400	19	36	0.28	0.58
C11 (GPa)	318 g	381	314	348	361	337
C12 (GPa)	103 g	109	99	102	112	101
C44 (GPa)	124 g	118	111	115	128	123
ɛ∞	5.35 h	6.34	6.34	6.30	5.08	5.20

^aLow temperature neutron measurements extrapolated to 0 K, reference [54]. ^bReference [117]. ^cSound velocity measurements at 298 K, reference [118]. ^dBrillouin scattering at 93 K, reference [119]. ^eReference [120]. ^fHigh-angle x-ray diffraction measurements, extrapolated to 0 K, reference [108]; back-corrected value is ${a}_{\text{ref}}^{\ast }=3.886$ Å. ^gSound velocity measurements. Values extracted at 273 K, reference [121]. ^hPermittivity measurements at room temperature, reference [122].

Meanwhile, the measured value a₂₉₈ = 3.905 Å provides a framework for our discussion of the functional-specific SrTiO₃ response characterizations (table 3). The elastic and dielectric constants are measured at room temperature, while we compute these properties at a₀. We shall, below, assign more trust in a given functional-specific response characterization when |a₀ − a₂₉₈| is small.

Table 3 reports a summary of the BO characterizations of BaZrO₃ and SrTiO₃ as computed using LDA, PBE, CX, HSE, and CX0P. For BaZrO₃ we find that the high-frequency (electronic component of the) dielectric constant _∞ is best described by the three non-hybrids, LDA, PBE, and CX, while the low-temperature data suggests that CX and CX0P provide the best characterizations for elastic constants. The results for _∞ (here obtained at the BO lattice constants) concur with characterizations reported in our previous BaZrO₃ study (that also tracked the impact of zero-point and thermal vibrational effects [55])—it is when adding the vibrational contributions in (0), that CX0P gains a performance edge for BaZrO₃ [54, 55]. For the elastic constants there are some scatter in the experimental data. However, if we rely on the low temperature measurements [119], we find that CX and CX0P are most accurate; we note that we must expect a lowering of these dielectric-constant results by the vibrational-driven expansion of the lattice.

For SrTiO₃ we find that CX0P provides the best BO characterization of the elastic constants overall but PBE and CX are also accurate. We repeat that these measurements are done at room temperature (where a₂₉₈ = 3.905 Å) while we compute them at the BO lattice constants. Taking the difference a₀ − a₂₉₈ into consideration we expect that the CX0P (PBE) values would soften (harden); the CX characterization can be directly compared with room-temperature measurements since they are provided for a cubic unit cell that corresponds almost exactly to the room-temperature structure (having a vanishing a₀ − a₂₉₈ difference); this fact increases the value of the CX accuracy that we have here documented. Meanwhile, we find that the characterizations of _∞ are similar among the hybrids and among the non-hybrid functionals (as in the BaZrO₃ case). For SrTiO₃ we also find that HSE and CX0P are significantly closer to the experimental results for _∞ [122].

Some of us have previously helped to establish (by combining experiment and theory) that BaZrO₃ remains cubic all the way down to T → 0 [54] and that CX0P stands out (over PBE/HSE and CX) in correctly describing measurements of the overall temperature variations in BaZrO₃ properties [54, 55]. Our additional results for the BaZrO₃ R₂₅ and Γ₁₅ modes (reported in table 3 at the BO structures), corroborate this conclusion. For BaZrO₃ we also find that use of PBE, HSE, and CX0P functionals directly leads to predictions of both ferroelectric stability (positive Γ₁₅ values) and AFD stability (positive R₂₅ values). LDA and CX characterizations do identify a potential BaZrO₃ AFD instability (imaginary R₂₅ values) at the BO structure. However, the magnitude of the R₂₅ mode value is so small in the CX characterization that zero-point dynamics and thermal effects stabilize BaZrO₃ at ambient conditions, as documented in reference [55].

We argue that the key to discussing functional success for SrTiO₃ is the answer to this question: 'Which of the functionals, if any, are consistent with the experimental SrTiO₃ finding of an actual low-temperature instability, driven by the SrTiO₃ AFD mode'? The answer should identify the functionals that are robust and can make the most relevant predictions. The answer helps us decide if we can trust the accuracy that is suggested in table 3.

The SrTiO₃ results for the R₂₅ and Γ₁₅ modes (in table 3) challenge us to pursue a detailed discussion of phase stability in materials. For the SrTiO₃ AFD R-mode and for the Γ modes, we find soft modes for all functionals (except in the LDA-Γ₁₅ description) in their description of the cubic phase. Our PBE and HSE results for the soft Γ₁₅ mode are in fair agreement with the values (PBE at i14 meV and HSE at i9 meV) reported in reference [124] and with the HSE value (i12.5 meV) value reported in reference [125]. We note that the BaZrO₃ Γ₁₅ mode is documented to be sensitive to the energy cutoff [55] and that the present SrTiO₃ results are obtained at considerably larger energy-cutoff values than in those previous studies.

The SrTiO₃ results for the R₂₅ mode identify incipient instabilities for all investigated functionals. That is, in all characterizations we find that the functionals have modes that could correspond to a possible structural transformation. The first question to address in our analysis is whether those predictions survive in the presence of the expected vibrationally-driven lattice expansion.

Figure 4 shows calculations of the potentially unstable R AFD (dashed curves) and Γ modes (solid curves) in cubic SrTiO₃ for different assumed lattice constants for the set of XC functionals LDA, PBE, CX, HSE, and CX0P. The dotted horizontal line identifies the zero-ω² value to delineate stability from a potential for instability; this is the measure of SrTiO₃ stability that applies in a phonon-level description, i.e., if we can ignore the stabilization that may emerge with a more complete account of the ZPE dynamics, see appendix B. The large squares show the position of the optimal BO lattice constant a₀, for each investigated functional; the vertical dotted line shows the experimental lattice constant as extrapolated to zero temperature [123]. For comparison, we note that the distance from the HSE (CX0P) BO lattice constant is about half (twice) the expected ZPE-lattice expansion 0.008 Å.

Zoom In Zoom Out Reset image size

We find, on the one hand, that HSE delivers a description that is close on the lattice constant, but on the other hand, yields an R-mode ${\omega }_{R}^{2}$ value that is barely negative at the optimal BO structure (let alone after the expected expansion). Quantum fluctuations are therefore expected to easily compensate and prevent the occurance of an AFD-driven phase transition in an HSE-based materials modeling of SrTiO₃ (and the expectation is substantiated by the discussion we present below for CX0P). Worse, the HSE Γ-mode ${\omega }_{{\Gamma}}^{2}$ value rapidly decreases beyond the BO lattice constant so that an HSE-based modeling seems to instead imply a possible ferroelectric SrTiO₃ behavior, which is again in conflict with experimental observations [107, 110].

Meanwhile, the use of CX0P gives ${\omega }_{R}^{2}$ values that are more negative than those for HSE at the native CX0P lattice constant. It also yields a modal description with an improved resistance toward a Γ mode (ferroelectric) instability, figure 4. However, CX0P is still in conflict with experimental observations: when the CX0P input (for ${\omega }_{R}^{2}$) is adjusted to the estimate for the T → 0 cubic lattice constant (vertical dashed line), it gives again only a very weak AFD-mode instability. Use of CX0P does not lead to the prediction that the AFD mode drives an actual low-temperature transformation in SrTiO₃.

To proceed with checks of general functional performance, we next compute the set of SrTiO₃ (and BaZrO₃) potential-energy surfaces (PES) for AFD-type deformations, for example, as shown for CX and CX0P in figure 5. In practice, we compute (for all chosen functionals) the energy cost associated with making a distortion in a 40-atom unit (top right panel of figure 3) reflecting a Glazer angle [126, 127] rotation around the z-axis, in steps of approximately 0.14 degrees.

Zoom In Zoom Out Reset image size

Landau-expansion theory [44] offers a natural framework for discussing DFT-based predictions of phase stability, for example in SrTiO₃. That is, the difference between finding a potential AFD instability (identified by having a $\bar{\kappa }< 0$ value) and predicting an actual instability (with a chosen functional) can be resolved by analyzing the nature of the deformation mode within a fourth-order Landau model. The model is defined by equation (B3) in appendix B and its use was illustrated for pressure-induced transitions BaZrO₃ in references [54, 55]. The Landau parameters, $\bar{\kappa }$ and $\bar{\alpha }$, are set from fitting to the underlying DFT results, for example, figure 5 for SrTiO₃. We obtain such descriptions for each perovskite and functional investigated, as also reported in table 3; the quality and the consistency of the Landau modeling (given by $\bar{\kappa }$ and $\bar{\alpha }$) is confirmed by checking against the R-mode description that arises in a harmonic approximation (as listed in table 3) using the approximation ${\omega }_{R}^{2}\approx \bar{\kappa }$. The Landau model is cast in a Hamiltonian form (appendix B) and we reveal the nature of the relevant deformation mode by solving the one-particle Schrödinger equation in the potential given by the PES, using the finite difference method.

The top (bottom) panel of figure 5 shows CX0P (CX) results for the total energy variation that occurs in SrTiO₃ when we track the AFD-type distortion while keeping the unit-cell lattice constants fixed at the optimal CX0P (CX) value (in line with the BO approximation). We note that large negative $\bar{\kappa }$ values of the Landau description correlate with deep double-well potentials. This is the case that emerges for the LDA and CX (and to some extent for the PBE) characterizations of SrTiO₃. However, the total-energy variation (and hence the potential for the effective AFD-model Hamiltonian, appendix B) is shallow for CX0P.

Importantly, the panels of figure 5 show the eigenlevels ω_i that emerge in our AFD-mode Landau modeling for SrTiO₃, as based on CX0P and CX input, see appendix B. The dashed horizontal lines depict the results for the eigenlevels as obtained under the assumption that the system is actually deformed; in this case the AFD dynamics occurs as trapped in one of the two displaced harmonic oscillators described by Q = ±Q₀ (as illustrated by dashed parabolae). The pair of solid horizontal lines—in each panel—show the eigenlevels for the AFD-modal dynamics as described under the assumption that there is no relevant dephasing of the modal double-well dynamics [46, 47, 49]. That is, this double-well eigenvalue description of the vibration is provided under the condition that the mode retains coherence and thus exists on both sides of the central barrier at Q ∼ 0.

To set up a simple stability criterion, we consider the SrTiO₃ (and general incipient-instability) case as an inelastic tunneling problem [46–50]. We interpret the 2Δ splitting of the eigenmodes for the AFD modal dynamics as arising from a hybridization characteristic for each of the functionals, as further discussed in appendix B. The computed magnitude of Δ reflects the rate of interwell tunneling and sets the scale of the characteristic dwell or tunneling time:

$$\begin{equation}{\tau }_{R}=\hslash /{\Delta}E.\end{equation} \tag{ 12 }$$

We compare these functional-specific results with an assumed inelastic-scattering or dephasing time τ_scat; a crude indicator for actually predicting a T → 0 phase transition is then

$$\begin{equation}{\tau }_{R}\gg {\tau }_{\text{scat}}.\end{equation} \tag{ 13 }$$

In essence, this criterion expresses the competition between the phase-coherence life time and the tunneling dynamics of the mode that could possibly drive a transformation. The criterion equation (13) is further discussed in appendix B. This appendix also motivates the use of τ_scat ∼ 1 ps as a natural delineation between the presence or absence of an actual instability in a given DFT-based modeling.

In table 3 we list the tunneling times τ_R as extracted for the set of functionals, using equation (12). We find that LDA is characterized by very long (ns) dwell or tunneling times, while CX and PBE give moderately long times (36 and 19 ps, respectively). Finally, the table reports that CX0P and HSE are characterized by short dwell times (less than 0.6 ps). With an assumed dephasing time on the order of 1 ps, our stability analysis suggests that the CX-based (CX0P-based) modeling predicts (does not predict) an actual distortion at T → 0.

In summary, we find that the CX provides the best overall description for SrTiO₃. Our conclusion is based on the observations that CX has a strong overall performance for properties that we can directly assert in DFT (table 3) and that the CX predictions for the AFD-mode behavior are consistent with experimental observations, unlike for CX0P. The finding of a CX performance edge over CX0P for modeling SrTiO₃ is in contrast to what we recently documented for BaZrO₃ [54, 55].

The GMTKN55 is a suite of benchmarks of broad molecular properties that also contain a range of DNA-relevant benchmark sets. For a first test of the CX, CX0P and AHCX performance on molecules it is relevant to consider the S66 set within the GMTKN55 suite [87]. The S66 is a set that broadly reflects noncovalent interactions and contains some nucleobase interactions.

More biorelated checks can be extracted by also computing MADs of our XC functional descriptions (relative to coupled-cluster CCSD(T) values [87]) for the PCONF21 set of peptide conformers, the amino20X4 set of amino-acid interaction energies, the UPU23 set of RNA backbone conformer energies, the SCONF set of sugar conformers, and the WATER27 set. We note that the sugar behavior also helps define the DNA back bone. Overall, we have thus extracted (from GMTKN55) a subset that focuses on biochemistry-related systems: nucleobases, amino acids, peptides, as well as RNA (and in part DNA) back bone properties.

Our testing setup is similar to that used in reference [15], but here we include also the WATER27 benchmarking set by computing the energy of the OH⁻ ion in a smaller 12 Å cubic cell. This allows us to circumvent adverse convergence impact of self-interaction errors in this negatively charged radical [15, 152].

Figure 7 shows a performance comparison for the CX-based consistent-vdW-DF tool chain, CX, CX0P, and AHCX. The performance of the consistent-vdW-DFs (CX/CX0P/AHCX) is strong overall on molecules [15] and very strong for most of the here-investigated bio-relevant problems. The performance is clearly better than, for example, that of the vdW-DF1 version [17].

Zoom In Zoom Out Reset image size

Table 4 summarizes our full comparison of biomolecule performance. The top three rows quantify the figure 7 performance overview for the CX-based tool chain, reporting MAD values (in kcal mol⁻¹) for the specific benchmark sets as well an average MAD measure (obtained by assuming equal weights among the 6 benchmark results for every functional). The table permits a quantitative comparison of CX, CX0P, and AHCX performance with benchmarking results obtained for other members of the vdW-DF family and for two dispersion-corrected functionals, revPBE + D3 [33, 85, 86] and HSE + D3 [12, 86]. The latter are closely related to the PBE-based semilocal tool chain and we note that revPBE + D3 is identified as the top performer of the Grimme-D3 corrected GGAs [87]. revPBE + D3 is a generalist that also does well for the here-considered biomolecular subset of the full GMTKN55 molecular testing suite [87]. The CX has a similar generalist status [15, 25], and table 4 shows that CX matches the revPBE + D3 performance also in this more targeted survey. Moreover, we find that the hybrid forms, CX0P and AHCX, improve the accuracy for biomolecular problems, as asserted in this test.

Interestingly, table 4 identifies the vdW-DF2 as the best overall performer, primarily because of its ability to accurate describe the energies involved in the WATER27 processes. From our broader molecular surveys, summarized in references [15, 25], we find that vdW-DF2 is only a fair competitor to the CX performance on noncovalent interactions, but it is exceptionally strong in a select few benchmarks, such as WATER27. This fact deserves a separate exploration, that we generally defer. However, for the ferroelectric polymer-crystal characterization (below), we include results for both vdW-DF2 and CX, comparing also with previous theory results [147, 149].

The WATER27 benchmark set constitutes a challenge for most XC functionals [25, 87]. Even here we find that the consistent-vdW-DF tool chain delivers a robust description, with a MAD value of 2.8 kcal mol⁻¹ for AHCX and almost as good for the nonhybrid CX form.

Overall, this survey suggests that CX and its tool chain are useful for determining interaction energies in biochemistry and, we expect, in both bio- and synthetic polymers.

Additionally, we test the CX performance and robustness against recent higher-level calculational results on DNA intercalation [128]. That study identifies a set of relevant frozen (reference-coordinate) geometries, figure 6, for which it also provides coupled cluster CCSD(T) reference energies of the energy gain by molecule insertion (ignoring elastic-energy costs). We use those results for an extra test of the CX performance, because the applied DNA modeling circumvents the need for a detailed study of the effects of counter ions and water: we can directly compare our DFT results against listed reference energies at specified coordinates [128].

The basic idea is to consider two models of a DNA base-pair segment, namely one where the backbone is protonated (effectively placing one extra electron per phosphor group on the back bone structure), and one where the back-bone is further eliminated. In reference [128] the atomic positions of the three intercalants, along with those of the immediate surrounding DNA structure, were extracted from the protein database (PDB) [153]. The structures used for the CCSD(T) results were truncated to the base pairs above and below the intercalant, plus the part of the sugar-phosphate backbone that connects them (model 'B', top panel of figure 6), or without the backbone (model 'A'). The interaction energies were calculated using the focal point approach [154] for extrapolated CCSD(T) results. There are reference energies (and structures) for both models with 3 intercalants, that are all effectively nearly flat, figure 6, for the parts that are inserted in the DNA; in addition, there are also reference energies for a variant process, where the intercalant is itself protonated.

Table 5 summarizes our CX results for the energy gains by DNA intercalation for the three intercalants and in model 'A' and 'B'. The comparison is made against the CCSD(T)-extrapolated results from reference [128], and with their B3LYP-D3, and HF-3c results. The three intercalants have PDB codes 1K9G, 1DL8 and 1Z3F and are here (and in reference [128]) denoted '1', '2', and '3' (see top and middle panels of figure 6). The latter molecule includes a nitrogen atom that can be protonated and this is also included in the set of calculations, the protonated molecule is denoted '3⁺'.

Inspecting the numbers in table 5, we see that CX yields results that are close to those from the CCSD(T) calculations, also for the protonated molecule ('3⁺'). This holds regardless if the DNA backbone is included (model 'B') or not (model 'A'). The largest deviation for CX is seen for molecule '1' (in both DNA models), with a 1.8–2.4 kcal mol⁻¹ difference from the CCSD(T)-energies. All other deviations are less than 0.7 kcal mol⁻¹, and MAD is 0.82 kcal mol⁻¹ for CX with respect to CCSD(T). This can be compared to the B3LYP-D3 deviations from CCSD(T), that show a slightly smaller MAD (0.54 kcal mol⁻¹) on the set of intercalate structures. However, B3LYP-D3 (and HF-3c) are hybrid calculations, and in plane-wave codes molecular hybrid calculations are found to be up to 30 times slower compared to regular functionals, like CX, for similar system sizes [15].

Turning to a comparison with results for the minimal-basis method HF-3c [128], we see that the MAD value for HF-3c, at 2.00 kcal mol⁻¹, is more than double that for CX. In other words, CX competes well with the best hybrid results of reference [128], and offers a path to acceleration that gives improved predictability compared to the HF-3c minimal-basis-set approach.

Finally, for our soft-matter application study, we characterize primarily the β crystalline form of the PVDF system, while also comparing with the so-called γ forms, figure 6. We predict the relaxed structure and ferroelectric response of perfect crystals of β- and γ-PVDF, while comparing with experiments and other theory results when possible. Both forms can be synthesized but, to the best of our knowledge, large single-crystal samples, with long-range order, do not yet exist. Theoretical predictions are thus motivated and we here seek to provide primarily a CX characterization in a 2 × 2 × 8 Monkhorst-Pack grid sampling of the Brillouin zone. We compute the single-chain energy E_chain in a unit-cell that has twice the lateral (non-chain) extensions than what holds in the experimental characterizations.

In analogy with equation (11), we compute the polymer-crystal cohesive energy,

$$\begin{equation}{E}_{\text{coh}}^{\text{pol}}(a,b,c)={E}_{\text{crystal}}-2{E}_{\text{chain}}\end{equation} \tag{ 14 }$$

for a series of unit-cell lattice constants in the β-PVDF form as well as for the motifs of the γ-PVDF form. We use the stress-vdW-DF implementation in QUANTUM ESPRESSO, although we did not have to here rely on the new spin-stress extension.

Figure 8 shows an overview of the structure search and potential-energy landscape for deformations of β-PVDF that we have computed using CX (as well as in vdW-DF1 and vdW-DF2). We first establish the energy dependence of the along-chain lattice constant (using constrained variable-cell calculations) as shown in the left panel. Then, at the optimal value c₀ we can extract the overall structure characterization, identified by the star shown in the right panel. In this panel, we furthermore report computed CX results for the energy variation ${E}_{\text{coh}}^{\text{pol}}({c}_{0},a,b)$ and a contour mapping of a fourth-order polynomial fit, tracking the energy variation in the soft interchain directions a and b [104, 139, 148].

Zoom In Zoom Out Reset image size

Table 6 summarizes the CX structure characterization along with those obtained using vdW-DF1 and vdW-DF2 and those found in literature. We note that the c₀-lattice constant is set by covalent interactions and that all functionals are overall in fair agreement of c₀. The functional characterizations differ primarily by the predictions of the soft lattice constants a and b. It is therefore natural to use the contour plot to illustrate the functional variations on PVDF structure determinations, as also done in the (right panel) of figure 8.

We find that CX results obtained with the above-described constrained-fit procedure align with those obtained in a constraint-free variable-cell optimization. However, that is true only when we start the system description of the polymer crystal close to the actual ground-state structure. In table 6 we give an extra decimal in the reporting to facilitate this comparison.

We also find that CX provides a highly accurate prediction of the optimal c lattice constant for the β-form relative to experiments; that is, the description of unit-cell extension along the polymer chains is almost spot-on the experimental observation. In the other directions, we find an overestimation of the a lattice constant but an underestimation of the b lattice constant. The PBE0 and vdW-DF1 results are in closer agreement with the b lattice constant, but there the c lattice constant is significantly overestimated.

The study of the γ-PVDF crystal offers additional opportunities for a theory-experiment comparison. These crystal forms are only nearly orthorhombic, characterized by a tilt angle between the a–b basis plane and the c axis that is (as for β-PVDF) aligned with the polymer strains. The comparison is complicated by the fact that there are two conformers, denoted γ_u and γ_d , see figure 6.

Table 7 shows the results of a vdW-DF1, vdW-DF2, and CX structure characterization, along with experimental observations for γ-PVDF [157]. These reference values are obtained at room temperature in systems that are known to contain a mixture of the two motifs, i.e., both γ_u and γ_d . Thus, one would expect the volume Ω, the tilt-angle, and the lengths of the set of unit-cell vectors to be a linear combination of γ_u and γ_d values. Moreover, polymers are known to exhibits a significant temperature expansion [150]. Consequently, a good description should have a predicted volume that is smaller than the volume measured at 300 K.

Table 7. Unit-cell parameters for γ-phases of PVDF. We compare present results as well as literature values. The experimental values rely on a sample that contains a mixture of the up and down configurations, characterized at T = 300 K. The unit cell is nearly orthorhombic, with a small tilt of the along-strain axis c and the a–b basis plane.

Functional	$a(\mathring{\text{A}})$	$b(\mathring{\text{A}})$	$c(\mathring{\text{A}})$	${\Omega}({\mathring{\text{A}}}^{3})$	∠ac(°)
γ d phase
PBE + D2 a	9.09	4.96	9.33	420.2	89.9
vdW-DF1 a	9.38	5.16	9.46	457.0	90.0
vdW-DF2 a	9.12	5.04	9.45	434.3	89.8
vdW-DF1	9.41	5.08	9.50	455.0	90.0
vdW-DF2	9.17	4.99	9.43	431.9	90.2
CX	9.31	5.02	9.60	449.3	90.0
γ u phase
PBE + D2 a	9.38	4.78	9.24	414.2	93.3
vdW-DF1 a	9.71	4.95	9.34	447.7	94.4
vdW-DF2 a	9.41	4.85	9.33	424.8	94.5
vdW-DF1	9.60	4.98	9.31	444.3	96.5
vdW-DF2	9.34	4.86	9.31	421.4	96.0
CX	9.36	4.84	9.32	421.8	94.6
Experiment b
γ u & γ d	9.67	4.96	9.20	440.65	93

^aReference [149]; dispersion correction 'D2' from reference [156]. ^bReference [157].

We find that the computed structures for the γ_d motif differ from the room-temperature experimental characterization, for all vdW-DF functionals. The tilt angle differs and so does the prediction for the along-chain extension (length of c). The CX characterization of the γ_d motif is overall not as close to the measurement as vdW-DF2. On the one hand, this is noteworthy because CX is in related systems found accurate on molecular-crystal bonding and structure, as exemplified for β-PVDF, table 6, and in references [56, 150]. On the other hand, it is not clear that the sample in the experiment contained a large component of γ_d motif.

Assuming instead that the structure of the polymer sample is dominated by the γ_u motif, the set of vdW-DF predictions are closer to the measured data. Here, CX describes a unit-cell tilt angle that is in good agreement with the measured data. Also, both vdW-DF2 and CX are now found to give a unit-cell description that is 5% smaller than the room-temperature measurements and the CX lattice constant is now in good (within 1%) agreement with the measured c-axis extension.

Overall, the CX is found accurate on structure for β-PVDF and consistent with the mixed-γ-motif measurements.

For predictions of the PVDF polarization, we follow the ideas presented by Johnsson and co-workers [147, 149] while using the QUANTUM ESPRESSO implementation of the modern theory of polarization [89–94]. That is, we use a modified cell with one PVDF chain rotated 180°, to establish nonpolar forms for both the γ and β phases and proceed with calculations that trace the effects of electron displacements in the bandstructure description of the polymer crystals [94]. For the γ_u/d phases, the nonpolar nature of such reference states is discussed in reference [149]; here we focus on describing the steps that are needed to accurately determine the (larger) spontaneous polarization that results in the β phase.

First, we determine the lateral position of the axis of along-chain rotation (for each of the two chains in the unit cell) by projecting the carbon atoms onto the lateral plane defined by the a and b lattice vectors. Second, we rotate one chain around its axis and confirm the nonpolar nature of the resulting reference structure [147]; specifically, we use QUANTUM ESPRESSO to compute the polarization that we find to be an integer times the so-called polarization quanta [94] (in this case e/ac) for the response along the b direction. Third, we track the evolution of the computed polarization as we incrementally rotate the chain back toward the actual PVDF β-phase structure, noting shifts among multiple branches of the polarization description in the Berry-phase formulation [89, 90, 93]. Finally, we extract the spontaneous polarization result, correcting for the impact of these polarization-branch jumps [94].

Table 8 summarises the outcome of these polarization-response surveys comparing also to previous vdW-DF1 and vdW-DF2 results [147, 149] for the γ and β phases. We note that our calculations are performed under the assumption of having a perfect crystal at the optimal structure (as computed for each of the vdW-DFs). We provide these theoretical characterizations not as a performance assessment but as an application example.

We find that CX characterizations are consistent with those provided with vdW-DF1 and vdW-DF2 for all three phases of PVDF. We also find that our vdW-DF1 and vdW-DF2 results in turn are in fair agreement (at least for the β phase) with previous theory characterizations [147, 149]. Comparison with experimental results (also listed) is difficult because of the challenges in securing fully crystalline samples, and polarization measurements will be affected by compensating responses arising in different sample regions. Nevertheless, we confirm that the β-PVDF form has the highest limit on the polarization response of the investigated phases.

In summary, we find that with a realization of a β-PVDF single crystal form, this ferroelectric polymer has a spontaneous polarization that is significantly larger than what is seen in present measurements [135]. The β phase may eventually serve as a good organic ferroelectrics.