Institut de Química Computacional i Catàlisi and Departament de Química,
Universitat de Girona, 17003 Girona, Spain

# RPA, an accurate and fast method for
the computation of non-linear optical properties

Pau Besalú-Sala  Department of Chemistry and Pharmaceutical Sciences,
Amsterdam Institute for Molecular and Life Sciences (AIMMS) and Vrije
Universiteit Amsterdam, De Boelelaan 1083, 1081 HV Amsterdam, The Netherlands
Fabien Bruneval  Université Paris-Saclay, CEA, Service de recherche en
Corrosion et Comportement des Matériaux, SRMP, 91191 Gif-sur-Yvette, France
Ángel José Pérez-Jiménez  Department of Physical Chemistry, University of
Alicante, E-03080 Alicante, Spain
Juan Carlos Sancho-García  Department of Physical Chemistry, University of
Alicante, E-03080 Alicante, Spain  Mauricio Rodríguez-Mayorga  Department of
Physical Chemistry, University of Alicante, E-03080 Alicante, Spain
<EMAIL_ADDRESS>
###### Abstract

The accurate computation of non-linear optical properties (NLOPs) in large
polymers requires accounting for electronic correlation effects with a
reasonable computational cost. The Random Phase Approximation (RPA) used in
the adiabatic connection fluctuation theorem is known to be a reliable and
cost-effective method to render electronic correlation effects when combined
with the density-fitting techniques and the integration over imaginary
frequencies. We explore the ability of the RPA energy expression to predict
NLOPs by evaluating RPA electronic energies in the presence of finite electric
fields to obtain (using the finite difference method) static polarizabilities
and hyper-polarizabilities. We show that RPA based on hybrid functional self-
consistent field calculations yields as accurate NLOPs as the best-tuned
double-hybrid functionals developed today, with the additional advantage that
RPA avoids any system-specific adjustment.

##  1 Introduction

The development of lasers in the 1960s led to the blossoming of non-linear
optics (NLO). Nowadays, many technologies are based on non-linear optical
effects triggered by the non-linear optical properties (NLOPs) of materials.
For example, NLOPs have applications in optical signal processing 1, ultra-
fast switches 2, sensors 3, and laser amplifiers 4, among others 5. From the
computational perspective, the prediction of static NLOPs for materials
requires the evaluation of (hyper)-polarizabilities produced with accurate
methods 6, 7, 8 able to describe the response of the many-body wavefunction
$\Psi$ and in particular of the electronic density, $n(\mathbf{r})=N\int
d\mathbf{r}_{2}\ldots
d\mathbf{r}_{N}\Psi^{*}(\mathbf{r},\mathbf{r}_{2},\ldots,\mathbf{r}_{N})\Psi(\mathbf{r},\mathbf{r}_{2},\ldots,\mathbf{r}_{N}),$
to the external fields produced by the incident light beams (see for example
Ref. 9). The classic definition of the molecular static electric
(hyper)-polarizabilities comes from the Taylor series expansion of the field-
dependent dipole moment 10, 11, 12, 13 ${\boldsymbol{\mu}}({\bf F})$ or energy
14 $E({\bf F})$. In the case of the $E({\bf F})$, its components read as

$E({\bf F})=E({\bf
0})-\sum_{i}\mu_{i}F_{i}-\frac{1}{2!}\sum_{ij}\alpha_{ij}F_{i}F_{j}-\frac{1}{3!}\sum_{ijk}\beta_{ijk}F_{i}F_{j}F_{k}-\frac{1}{4!}\sum_{ijkl}\gamma_{ijkl}F_{i}F_{j}F_{k}F_{l}+\dots$
(1)

with $i,j,k$ and $l$ being any Cartesian component (i.e. $x,y$ or $z$) and the
expansion coefficients,

$\displaystyle\mu_{i}=-\frac{\partial E}{\partial F_{i}}\Big{|}_{{\bf F}={\bf
0}}\;,\;\alpha_{ij}$ $\displaystyle=-\frac{\partial^{2}E}{\partial
F_{i}\partial F_{j}}\Big{|}_{{\bf F}={\bf
0}}\;,\;\beta_{ijk}=-\frac{\partial^{3}E}{\partial F_{i}\partial F_{j}\partial
F_{k}}\Big{|}_{{\bf F}={\bf 0}}\;,$
$\displaystyle\;\textrm{and}\;\;\gamma_{ijkl}$
$\displaystyle=-\frac{\partial^{4}E}{\partial F_{i}\partial F_{j}\partial
F_{k}\partial F_{l}}\Big{|}_{{\bf F}={\bf 0}},$ (2)

being the dipole moment, the static polarizability, the first hyper-
polarizability, and the second hyper-polarizability, respectively. The
(hyper)-polarizabilities can be alternatively obtained directly from the
derivatives of the electronic density (see Ref. 15 and references therein).

Within density functional theory (DFT), the electronic energy is a functional
of the electronic density $E=E\left[n(\mathbf{r})\right]$ 16. In practical
applications using DFT, it is common to employ the Kohn-Sham (KS) scheme 17,
which is based on expressing the electronic many-electron wavefunction as a
single-determinant $\Phi$ and the density as
$n(\mathbf{r})=\sum^{occ}_{i}|\phi_{i}({\bf r})|^{2}$ where $\\{\phi_{i}\\}$
are the one-electron wavefunctions that solve the Kohn-Sham Hamiltonian (i.e.
$\widehat{H}^{KS}\phi_{i}=\varepsilon_{i}\phi_{i}$). Unfortunately, within KS
DFT the unknown exchange-correlation energy functional ($E_{xc}$) needs to be
approximated. Numerous $E_{xc}$ functional approximations proposed in the
literature 18 also depend explicitly on the $\\{\phi_{i}\\}$ and
$\\{\varepsilon_{i}\\}$, like the so-called (double)-hybrid functional
approximations. In consequence, the total energy for the (double)-hybrid
functional approximations including the contributions from the nuclei within
the Born-Oppenheimer approximation 19 reads

$\displaystyle E\left[\phi,\varepsilon\right]$
$\displaystyle=-\frac{1}{2}\sum^{occ}_{i}\int d\mathbf{r}\;\phi^{*}_{i}({\bf
r})\nabla^{2}_{\mathbf{r}}\phi_{i}({\bf r})+\int
d\mathbf{r}\;v_{\textrm{ext}}(\mathbf{r})n(\mathbf{r})+\frac{1}{2}\int
d\mathbf{r}\;d\mathbf{r}^{\prime}\;\frac{n(\mathbf{r})n(\mathbf{r}^{\prime})}{|\mathbf{r}-\mathbf{r}^{\prime}|}+E_{nn-
nf}$
$\displaystyle+aE_{x}^{\textrm{EXX}}\left[\\{\phi\\}\right]+(1-a)E_{x}^{\textrm{DFT}}\left[n,\nabla
n\right]+bE_{c}^{\textrm{PT}}\left[\\{\phi,\varepsilon\\}\right]+(1-b)E_{c}^{\textrm{DFT}}\left[n,\nabla
n\right],$ (3)

with $a$ and $b$ being the hybridization parameters, $E_{nn-nf}$ being the
nucleus-nucleus interaction plus any nucleus-external field interaction,
$v_{\textrm{ext}}(\mathbf{r})$ being the time-independent external potential
produced by the (bare) nuclei plus any electric field applied ${\bf F}$.
$E_{x}^{\textrm{DFT}}$ ($E_{c}^{\textrm{DFT}}$) is a DFT functional
approximation that accounts for the exchange (resp. correlation) often
depending also on the gradient of the electronic density $\nabla n$ 20,
$E_{x}^{\textrm{EXX}}$ is the exact exchange obtained using $\Phi$ (i.e. the
occupied orbitals), and $E_{c}^{\textrm{PT}}$ is a term approximated from
perturbation theory 21 (i.e. $E_{c}^{\textrm{PT}}\sim E_{c}^{\textrm{MP2}}$)
that accounts for electronic correlation effects With MP2 being the second-
order perturbation theory correction as proposer by Møller and Plesset
moller_pr1934. Unfortunately, in general, KS DFT functional approximations
fail to predict accurate NLOPs due to the delocalization error 23, the self-
interaction error 24, and the incorrect decay of the exchange-correlation
kernels 25. By partially correcting these issues, the so-called range-
separated functionals 26, 27, 28 such as LC-BLYP 29, CAM-B3LYP 30, among
others have shown to improve in the description of NLOPs w.r.t. to their non-
range separated counterparts. Further optimal tuning (OT) of range-separated
functionals can yield better NLOP predictions. Among all the (OT)-range-
separated functionals, the T$\alpha$-LC-BLYP functional approximation,
proposed by one of us, 15,Built using the most popular OT technique that
consists of enforcing the analogous of Koopmans’ theorem for DFT
koopmans1934zuordnung, janak_prb1978. has proven to be the best approximation
among the tested functionals for computing NLOPs for a set of different
organic and inorganic molecules. Stepping up to the highest rung that is
currently in use in Jacob’s ladder picture, 32, 33 we find the KS DFT
functional approximations with a certain amount of MP2 correlation energy
(a.k.a. double-hybrid functionals in the literature 34, 35, 36) such as B2PLYP
34, PBE0-DH 35, and PBE-QIDH 36, which have been proposed in an attempt to
improve the non-local character of the correlation energy at the expense of
increasing the computational cost. Let us stress that the ability of the most
recently developed double-hybrid approximations to predict NLOPs was still an
open question as of today.

Another method to produce electronic energies consists in using the random
phase approximation (RPA); setting $E_{c}^{\textrm{PT}}\sim
E_{c}^{\textrm{RPA}}$ in Eq. (1) with the $E_{c}^{\textrm{RPA}}$ given by the
adiabatic connection fluctuation-dissipation theorem expression, which can be
written using imaginary frequencies as 37, 38, 39

$E_{c}^{\textrm{RPA}}[\chi_{0}]=\frac{1}{4\pi}\int_{-\infty}^{+\infty}d\omega\,\mathrm{Tr}\left\\{v\chi_{0}(\mathrm{i}\omega)+\mathrm{ln}\left[1-v\chi_{0}(\mathrm{i}\omega)\right]\right\\}.$
(4)

where Tr indicates the trace, $v$ is the Coulomb interaction, and $\chi_{0}$
is the non-interacting (dynamic) polarizability that can be built in terms of
the KS DFT $\\{\phi_{i}\\}$ and $\\{\varepsilon_{i}\\}$. RPA has attracted
much attention from the scientific community 40, 41, 42, 43, 44, 45, 46, 47,
48 due to its ability to account for the so-called non-dynamic electronic
correlation effects 49 with a low computational cost (i.e. it scales as
$M^{3}$ with $M$ being the size of the basis set 50, 51, 52, 53). The usual
procedure for computing RPA electronic energies involves performing a KS DFT
calculation with a DFT functional approximation That may include a certain
amount of exact exchange. and setting $b=0$ in Eq. (1). Once the self-
consistent procedure is completed, the KS DFT $\\{\phi_{i}\\}$ and
$\\{\varepsilon_{i}\\}$ are inserted in Eq. (1) with $a=1$ and $b=1$ (i.e.
canceling all KS DFT contributions) to produce electronic energies. The good
ratio of accuracy and computational cost of RPA makes it a very appealing
approximation for the computation of NLOPs for large systems: it is the main
focus of this work.

##  2 Methodology

We have implemented in MOLGW code 55, 56 the homogeneous electric fields ${\bf
F}$, needed for the computation of the numerical (hyper)-polarizabilities in
Eq. (1). In this work, we have used finite differences computed with an in-
house program that utilizes the Rutishauser–Romberg scheme 57, 58 to reduce
the numerical errors (see Refs. 8 and 15 for more details). Furthermore, the
KS DFT double-hybrid functional approximations tested in this work (i.e.
B2PLYP, PBE0-DH, and PBE-QIDH) have been introduced in MOLGW code The latest
release of MOLGW code (version 3.2) incorporates these new features (download
from https://github.com/bruneval/molgw/releases/tag/v3.2).. The aug-cc-pVDZ
basis set 60, 61, 62 has been chosen to facilitate the comparison with
previous studies 15. In addition, the density-fitting technique has been
employed with the aug-cc-pVDZ-RI auxiliary basis set 63, 64 in all
calculations. The RPA correlation energy given in Eq. (4) has been computed
using the density fitting technique and a Gauss–Legendre quadrature with 40
imaginary frequencies, which has proven to ensure high accuracy in the
integration with errors lower than $10^{-6}$ a.u. The set of systems employed
in this work has been taken from Ref. 15, where coupled-cluster singles and
doubles with perturbative triples correction $[$CCSD(T)$]$ values were used as
reference. In particular, 50 representative molecular systems formed by atoms
of the second period and/or hydrogen (see Fig. 1) were chosen from Ref. 15,
including representative oligomers among the most typical $\pi$-conjugated NLO
polymers, such as all-trans polyacetylene (PA), all-trans poly-methyneimine
(PMI), and polydiacetylene (PDA). Lastly, some small organic and inorganic
molecules, and hydrogen chains (known to be challenging systems for the
calculation of NLOPs) are also included. Interestingly, the studied set had
been used for the construction of a long-range corrected range-separated
hybrid, called T$\alpha$-LC-BLYP, that has been OT to compute the NLOPs of the
molecules of this set, allowing us to conduct a strict assessment of the RPA
results.

Figure 1: Chemical structures of the molecular set employed in this study
including the orientation of the applied field. Subset A contains molecules
that present their odd energy derivatives equal to 0 (i.e.
$\mu_{z}=\beta_{zz}=0$) because they have inversion symmetry, while subset B
is formed by molecules with all their energy derivatives different from 0. The
geometries were taken from the supporting information of Ref. 15 (and are also
available in Ref. 65).

##  3 Results

The aim of this work is to scrutinize the performance of the RPA energy
functional given in Eqs. (1) and (4) to predict polarizabilities and hyper-
polarizabilities and decide whether RPA can be used routinely in the
computation of NLOPs thanks to its (relatively) low computational cost 50, 51,
52, 53. Also, the double-hybrid functional approximations (i.e. B2PLYP,
PBE0-DH, and PBE-QIDH) are considered to be the highest rung in Jacob’s ladder
picture; thus, it is time to appraise the performance of the most recently
developed approximations for predicting NLOPs. Hence, we have computed the
$\mu_{z}$, $\alpha_{zz}$, $\beta_{zzz}$, and $\gamma_{zzzz}$ values for the
systems introduced in the previous section using RPA, B2PLYP, PBE0-DH, and
PBE-QIDH. In general, our results show that the absolute errors are dominated
by $\gamma_{zzzz}$, where they present the largest deviations w.r.t. the
CCSD(T) reference values. 15 It is worth mentioning that the largest absolute
errors for each property were obtained for the large systems like PDA4 (see
the supporting information for more details) in line with the previous work
15.

Next, let us remark that RPA is computed in a non-self-consistent manner;
hence, we must assess the role of the amount of $E_{x}^{\textrm{EXX}}$ present
at the underlying self-consistent KS DFT level when using RPA. In previous
works 66, 67, it has been shown that large values of $a$ in Eq. (1) tend to
produce electronic energies that depend less on the starting point. Hence, we
also aim to confirm if large $a$ values in Eq. (1) employed at the KS DFT
starting point produce more accurate properties (i.e. NLOPs). For this end, we
have taken as reference the PBEh($a$) hybrid functional 68 and varied the
amount of $E_{x}^{\textrm{EXX}}$ (specifying this amount in parenthesis).

In Fig. 2, we have collected MAX% and MEAN% Notice that for RPA@PBEh(0.0), the
MAX% for $\alpha_{zz}$ is 38.47%, the MEAN% of $\gamma_{zzzz}$ is 96.85%, and
the MAX% of $\gamma_{zzzz}$ is 898.91%; hence, these values lie out of the
scale used. to be used as a yardstick for the predicted dipole moments and
(hyper)-polarizabilities using RPA@PBEh($a$) for different $a$ values. Our
results indicate that $\mu_{z}$ is better described using lower values of the
$a$ hybridization coefficient at the KS DFT level. Nevertheless, $\mu_{z}$ is
not strongly dependent on the amount of $E_{x}^{\textrm{EXX}}$, as it is shown
in Fig. 2. For $\alpha_{zz}$, the role of the amount of $E_{x}^{\textrm{EXX}}$
employed at the KS DFT level becomes more relevant. Clearly, a small amount of
exact exchange (such as $a=0.25$) employed at the KS DFT level for producing
$\\{\phi_{i}\\}$ and $\\{\varepsilon_{i}\\}$ is able to drop the relative
errors dramatically to approximately half of the values obtained when a non-
hybrid functional is employed as starting point (recalling that PBEh($0$) is
the pristine PBE functional). In addition, notice that the $a=0.25$ value
already produces MAX% and MEAN% that compare well with the same values
obtained with larger amounts of $E_{x}^{\textrm{EXX}}$ ($a>0.25$); therefore,
the role of the hybridization coefficient $a$ is reduced for $\alpha_{zz}$ as
soon as some exact exchange is already present at the KS DFT level. Moving to
the next derivative, $\beta_{zzz}$, it is remarkable that using $a=0.0$ (i.e.
the pristine PBE functional) or $a=0.5$ as starting points yields similar MAX%
values. Although, we observe an important decrease (of approximately 50%) in
the MEAN% using $a=1.0$, which indicates that large values of $a$ should be
preferred for predicting $\beta_{zzz}$ accurately. In sharp contrast, the
relative errors obtained for $\gamma_{zzzz}$ indicate that the role of the
amount of $E_{x}^{\textrm{EXX}}$ is crucial for its accurate description,
where a small amount of exchange is fundamental to drop the errors by about
80% w.r.t. the RPA@PBEh($0$) (i.e. RPA@PBE) ones. Our results show that an
optimal value of $0.5<a<1.0$ can be chosen for dropping the relative errors
when computing $\gamma_{zzzz}$ because the RPA@PBEh(0.5) and RPA@PBEh(1.0)
errors are larger than the RPA@PBEh(0.75) ones. In general, our analysis of
the dependence on the starting point reveals that some amount of
$E_{x}^{\textrm{EXX}}$ must be present at the KS DFT level to produce
acceptable $\\{\phi_{i}\\}$ and $\\{\varepsilon_{i}\\}$ to enter the RPA
energy expression. This is in line with the fact that DFT functionals with a
large amount of $E_{x}^{\textrm{EXX}}$ tend to be more accurate for NLOP
computations. The dependence on the starting point seems to be less evident
than for electronic energies reported in our previous works 66, 67 because for
each NLOP an $a$ value can be selected leading to accurate values for that
particular property, but not necessarily improving the description of the
other ones, with $\mu_{z}$ and $\beta_{zzz}$ being less dependent on the
starting point than $\alpha_{zz}$ and $\gamma_{zzzz}$. Lastly, we observe that
large values of $a$ (i.e. $a>0.75$) seem to be able to lead to a compromise
situation where $\mu_{z}$, $\alpha_{zz}$, $\beta_{zzz}$, and $\gamma_{zzzz}$
can all be reasonably well described with the same fixed $a$ value.

Figure 2: Mean relative errors (MEAN%) and maximum relative errors (MAX%)
obtained with RPA@PBEh($a$) for the predicted (hyper)-polarizabilities w.r.t.
CCSD(T) for different values of the $a$ hybridization coefficient employed at
the KS DFT level. Note: In the case of RPA@PBEh(0.0), the MAX% for
$\alpha_{zz}$ is 38.47%, the MEAN% of $\gamma_{zzzz}$ is 96.85%, and the MAX%
of $\gamma_{zzzz}$ is 898.91%. Thus, these values lie out of the scale.

From our previous analysis, we recognized that $a\in\left(0.5,1.0\right)$
could lead to an accurate prediction of $\beta_{zzz}$ and $\gamma_{zzzz}$ at
the same time, the two properties that usually DFT tends to fail in larger
proportion to evaluate. Indeed, we have found that a $a=0.85$ value leads to a
good compromise situation for both properties, which allows us to introduce
the RPA@PBEh(0.85) approximation for the description of NLOPs. In Table 1, we
have collected the mean absolute errors (MAE) in a.u., root-mean-square errors
(RMSE) in a.u., maximum error (ME) in a.u., MEAN%, and MAX% values obtained
with this method. As expected MAE, RMSE, and ME are dominated by large systems
(e.g. PDA4) and their error increases rapidly with the order of the
(numerical) energy derivative (i.e. they increase from $\mu_{z}$ to
$\gamma_{zzzz}$). From Table 1, we observe that for the MAE, RMSE, and ME that
the RPA@PBEh(0.85) approximation performs as well as the T$\alpha$-LC-BLYP
functional for $\mu_{z}$, where the largest deviations are obtained for the
PMI6 system. It is the best method for predicting $\alpha_{zz}$ with the best
description of PDA systems than any other method and with the largest
deviations obtained for PA systems. It deteriorates some for computing
$\beta_{zzz}$, where the largest deviation is that of the PMI6 system. And,
finally, it performs similarly to LC-BLYP functional for $\gamma_{zzzz}$ with
the largest errors obtained for the PDA systems. Nevertheless, we must stress
once more that for computing NLOPs of molecules of different sizes, the most
important errors are the relative errors (represented by MEAN% and MAX% in
Table 1), where the RPA@PBEh(0.85) approximation performs better than any non-
OT-range-separated functional (e.g. it is better than LC-BLYP) for dipole
moments presenting the largest deviations for small systems (i.e. HNO and CO).
Moving on to the computation of $\alpha_{zz}$, RPA@PBEh(0.85) is the best
method among all the tested ones, including OT-range-separated hybrids. For
all the functionals tested, the largest errors on $\alpha_{zz}$ were produced
by PA and PDA systems (in agreement with the results obtained for the MAE,
RMSE, and ME). Remarkably, despite PA/PDA causing the highest errors for
RPA@PBEh(0.85), those errors are lower compared to other functionals,
explaining the success of RPA@PBEh(0.85) on the polarizability calculations.
For $\beta_{zzz}$, RPA@PBEh(0.85) can also be considered the best method
because it shows the lowest maximum percentage error, but on the contrary to
$\alpha_{zz}$, in this case, there is not a specific family of molecules
causing the error but the overall. Finally, for $\gamma_{zzzz}$,
RPA@PBEh(0.85) is the second-best method behind the (system adapted)
T$\alpha$-LC-BLYP functional. This is expected as T$\alpha$-LC-BLYP was
specifically designed to compute $\gamma_{zzzz}$ on this particular set of
molecules 15. Interestingly, the largest deviations of RPA@PBEh(0.85) are
obtained for small systems (i.e. CO and H2O), while relative errors lower than
32% for all the polymers studied were obtained, and with excellent performance
for Hydrogen chains where relative errors are lower than 8%, which praises its
good performance for predicting $\gamma_{zzzz}$. It is worth stressing again
that RPA@PBEh(0.85) does not involve any extra re-optimization of parameters
for each system, which makes it a very competitive approximation for computing
NLOPs when compared with the best KS DFT functional currently in use 15.

Table 1: Mean absolute errors (MAE) in a.u., root-mean-square errors (RMSE) in
a.u., maximum error (ME) in a.u., mean relative errors (MEAN%), and maximum
relative errors (MAX%) for the predicted (hyper)-polarizabilities w.r.t.
CCSD(T) values._a_

Property |  | T$\alpha$-LC-BLYP _a_ | LC-BLYP | PBE-QIDH | RPA@PBEh(0.85)
---|---|---|---|---|---
$\mu_{z}$ | MAE (a.u.) | 0.12 | 0.16 | 0.12 | 0.11
RMSE (a.u.) | 0.20 | 0.25 | 0.20 | 0.19
ME (a.u.) | 0.53 | 0.67 | 0.52 | 0.50
MEAN% | 5.03 | 29.89 | 13.88 | 13.67
| MAX% | 67.11 | 174.32 | 46.36 | 52.61
$\alpha_{zz}$ | MAE (a.u.) | 9.61 | 29.70 | 12.32 | 5.34
RMSE (a.u.) | 17.54 | 99.63 | 27.64 | 10.89
ME (a.u.) | 543.23 | 548.02 | 199.39 | 192.99
MEAN% | 5.66 | 10.23 | 5.45 | 3.07
| MAX% | 13.15 | 82.29 | 20.38 | 9.73
$\beta_{zzz}$ | MAE (a.u.) | 101.08 | 52.30 | 96.35 | 142.10
RMSE (a.u.) | 276.54 | 382.42 | 140.09 | 348.27
ME (a.u.) | 1.11$\times$103 | 1.50$\times$103 | 5.55$\times$102 | 1.24$\times$103
MEAN% | 73.59 | 61.14 | 39.20 | 49.34
| MAX% | 433.48 | 413.13 | 258.14 | 121.06
$\gamma_{zzzz}$ | MAE (a.u.) | 1.88$\times$104 | 4.06$\times$104 | 1.51$\times$105 | 6.58$\times$104
RMSE (a.u.) | 5.85$\times$104 | 9.22$\times$104 | 6.35$\times$105 | 2.60$\times$105
ME (a.u.) | 3.22$\times$105 | 3.43$\times$105 | 4.16$\times$106 | 1.64$\times$106
MEAN% | 5.92 | 19.44 | 19.13 | 13.90
| MAX% | 29.32 | 50.73 | 78.45 | 49.50

_a_ CCSD(T) and T$\alpha$-LC-BLYP values taken from Ref. 15.

Among the double-hybrid functional approximations, the MAE, RMSE, and ME are
again dominated by large systems (see the supporting information for more
details). Nevertheless, focusing on relative errors we observe from Fig. 3
that PBE-QIDH outperforms B2PLYP and PBE0-DH for all polarizabilities and
hyper-polarizabilities. Actually, only for predicting dipole moments the
PBE0-DH functional is the best one (closely followed by PBE-QIDH).
Interestingly, the PBE-QIDH functional presents the largest amount of
$E_{x}^{\textrm{EXX}}$, which confirms that the role of the exact exchange is
fundamental for an accurate description of NLOPs. Since PBE-QIDH is the best
double-hybrid functional approximation employed in this work, we have
collected in Table 1 the MAE, RMSE, ME, MEAN%, and the MAX% obtained using
this functional approximation. Clearly, focusing on relative errors, PBE-QIDH
performs as well as RPA@PBEh(0.85) for $\mu_{z}$, it is very competitive for
computing $\alpha_{zz}$ with errors that lie close to the T$\alpha$-LC-BLYP
ones, it is a really good method for computing $\beta_{zzz}$ because its MEAN%
is the lowest one (39.2%), and, finally, it compares well with the best (non-
system adapted) range-separated functional (i.e. LC-BLYP) for predicting
$\gamma_{zzzz}$. In general, PBE-QIDH results indicate that it is also a
competitive method for computing NLOPs, but involving an increase in the
computational cost due to the evaluation of the MP2 correlation energy (that
can rapidly increase with system size). Finally, let us mention that the
increase in the complexity of the functional approximation due to the
inclusion of a double hybrid scheme needs an adequate selection of the
hybridization coefficients. In fact, the truly non-empirical design of PBE-
QIDH (without any fitted parameter), as defined by some of us 33, 36, leads to
a large amount of exact exchange present in this functional ($a=3^{-{1/3}}\sim
0.69336$). Consequently, the large amount of $E_{x}^{\textrm{EXX}}$ reduces
fundamental errors (such as the delocalization error, the self-interaction
error, among others) that can play an important role in the description of
properties using KS DFT approximations as we show in this work.

Figure 3: Mean relative errors (MEAN%) and maximum relative errors (MAX%) for
the predicted (hyper)-polarizabilities w.r.t. CCSD(T) for B2PLYP, PBE0-DH, and
PBE-QIDH.

To rationalize the results, let us recall that the functional expression used
in this work for the RPA electronic energy employs (as it is usually done)
100% of the exact exchange (i.e. $a=1$ in Eq. (1)); therefore, the well-known
errors 23, 15 that are responsible for the usual failures of KS DFT in the
description of NLOPs are drastically reduced. Moreover, the RPA is known to be
equivalent to the direct ring-coupled-cluster doubles method 70, which ensures
a reasonable description of the so-called dynamic electronic correlation
effects; approaching the coupled-cluster description for these effects. On the
other hand, the results obtained with the double-hybrid functional
approximations indicate that PBE-QIDH is the best double-hybrid KS DFT
approximation for computing NLOPs (among the ones tested in this work). The
good performance of PBE-QIDH can be justified due to a large amount of exact
exchange present in this functional ($a=3^{-1/3}\sim 0.69336$), while B2PLYP
contains a lower amount of $E_{x}^{\textrm{EXX}}$ ($a=0.53$) and PBE0-DH even
a lower value ($a=0.5$). Indeed, the increase in the amount of exact exchange
makes the PBE-QIDH reduce the delocalization error and the self-interaction
error, which leads to a very competitive performance in the computation of
polarizabilities and hyper-polarizabilities using this functional
approximation. Finally, comparing PBE-QIDH performance for predicting NLOPs
with that of RPA@PBEh(0.85), see Table 1, the latter overtakes the former with
less computational expense, which makes RPA@PBEh(0.85) a very attractive
approximation for computing these properties.

In summary, the RPA@PBEh(0.85) approximation is a simple, accurate, and fast
method that can be used routinely in the computation of non-linear optical
properties. Its low computational cost 50, 51, 52, 53 makes it a very
competitive approximation outperforming KS DFT functionals approximations. We
have also shown that PBE-QIDH functional performs well and is the best
(global) double-hybrid KS DFT functional approximation (currently available)
for computing NLOPs, at the expense of increasing the computational cost.
Since the PBE-QIDH functional approximation is a functional that does not
involve a system-dependent optimization, its universality also facilitates its
usage on different systems.

A. J. P.-J., J. C. S.-G., and M.R.-M. acknowledge the Ministerio de Ciencia e
Innovación de España for Grant No. PID2019-106114GB-I00. P. B.-S. acknowledges
the financial support received from the Vrije Universiteit Amsterdam. F.B. and
M.R.-M. acknowledge the financial support provided by the Cross-Disciplinary
Program on Numerical Simulation of the French Alternative Energies and Atomic
Energy Commission (CEA) (ABIDM project). We also acknowledge the computational
resources of the Institut de Química Computacional i Catàlisi of the
Universitat de Girona. Finally, the authors thank É. Brémond for the
discussions related to the implementation of double-hybrid KS DFT functionals.

Dipole moments, polarizabilities, and hyper-polarizabilities are given in the
supporting information (as a spreadsheet).

## Graphical abstract

Figure 4: TOC

## References

  * Astill 1991 Astill, A. Material figures of merit for non-linear optics. _Thin Solid Films_ 1991, _204_ , 1–17
  * Quémard et al. 2001 Quémard, C.; Smektala, F.; Couderc, V.; Barthelemy, A.; Lucas, J. Chalcogenide glasses with high non linear optical properties for telecommunications. _Journal of Physics and Chemistry of Solids_ 2001, _62_ , 1435–1440
  * Gounden et al. 2020 Gounden, D.; Nombona, N.; Van Zyl, W. E. Recent advances in phthalocyanines for chemical sensor, non-linear optics (NLO) and energy storage applications. _Coordination Chemistry Reviews_ 2020, _420_ , 213359
  * Mace and Adams 1987 Mace, D.; Adams, M. Aspects of Non-linear Optics. _Europhysics News_ 1987, _18_ , 48–51
  * Castet et al. 2022 Castet, F.; Tonnelé, C.; Muccioli, L.; Champagne, B. Predicting the Second-Order Nonlinear Optical Responses of Organic Materials: The Role of Dynamics. _Accounts of Chemical Research_ 2022, _55_ , 3716–3726
  * de Wergifosse et al. 2014 de Wergifosse, M.; Liégeois, V.; Champagne, B. Evaluation of the molecular static and dynamic first hyperpolarizabilities. _International Journal of Quantum Chemistry_ 2014, _114_ , 900–910
  * Champagne et al. 2018 Champagne, B.; Beaujean, P.; De Wergifosse, M.; Cardenuto, M. H.; Liégeois, V.; Castet, F. Quantum chemical methods for predicting and interpreting second-order nonlinear optical properties: From small to extended $\pi$-conjugated molecules. _Frontiers of Quantum Chemistry_ 2018, 117–138
  * Naim et al. 2023 Naim, C.; Besalú-Sala, P.; Zaleśny, R.; Luis, J. M.; Castet, F.; Matito, E. Are Accelerated and Enhanced Wave Function Methods Accurate to Compute Static Linear and Nonlinear Optical Properties? _Journal of Chemical Theory and Computation_ 2023, _19_ , 1753–1764
  * Zhenyu and Xin 2023 Zhenyu, Z.; Xin, X. Focal-Point Analysis to Achieve Accurate CCSD(T) Data Set References for Static Polarizabilities. _Journal of chemical theory and computation_ 2023,
  * Chopra et al. 1989 Chopra, P.; Carlacci, L.; King, H. F.; Prasad, P. N. Ab initio calculations of polarizabilities and second hyperpolarizabilities in organic molecules with extended. pi.-electron conjugation. _J. Phys. Chem._ 1989, _93_ , 7120–7130
  * Buckingham and Orr 1967 Buckingham, A.; Orr, B. Quarterly Reviews. Chem. Soc. London. 1967; p 195
  * Orr and Ward 1971 Orr, B.; Ward, J. Perturbation theory of the non-linear optical polarization of an isolated system. _Molecular Physics_ 1971, _20_ , 513–526
  * Ditchfield et al. 1970 Ditchfield, R.; Ostlund, N.; Murrell, J.; Turpin, M. Comparison of the sum-over-states and finite perturbation theories of electrical polarizability and nuclear spin-spin coupling. _Molecular Physics_ 1970, _18_ , 433–440
  * Nakano et al. 1995 Nakano, M.; Shigemoto, I.; Yamada, S.; Yamaguchi, K. Size-consistent approach and density analysis of hyperpolarizability: Second hyperpolarizabilities of polymeric systems with and without defects. _The Journal of Chemical Physics_ 1995, _103_ , 4175–4191
  * Besalú-Sala et al. 2020 Besalú-Sala, P.; Sitkiewicz, S. P.; Salvador, P.; Matito, E.; Luis, J. M. A new tuned range-separated density functional for the accurate calculation of second hyperpolarizabilities. _Physical Chemistry Chemical Physics_ 2020, _22_ , 11871–11880
  * Hohenberg and Kohn 1964 Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. _Phys. Rev._ 1964, _136_ , B864–B871
  * Kohn and Sham 1965 Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. _Phys. Rev._ 1965, _140_ , A1133–A1138
  * Burke and Wagner 2013 Burke, K.; Wagner, L. O. DFT in a nutshell. _International Journal of Quantum Chemistry_ 2013, _113_ , 96–101
  * Combes et al. 1981 Combes, J.-M.; Duclos, P.; Seiler, R. The born-oppenheimer approximation. _Rigorous atomic and molecular physics_ 1981, 185–213
  * Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. _Phys. Rev. Lett._ 1996, _77_ , 3865–3868
  * Görling 1996 Görling, A. Exact treatment of exchange in Kohn-Sham band-structure schemes. _Phys. Rev. B_ 1996, _53_ , 7024–7029
  * 22 With MP2 being the second-order perturbation theory correction as proposer by Møller and Plesset moller_pr1934
  * Bryenton et al. 2023 Bryenton, K. R.; Adeleke, A. A.; Dale, S. G.; Johnson, E. R. Delocalization error: The greatest outstanding challenge in density-functional theory. _Wiley Interdisciplinary Reviews: Computational Molecular Science_ 2023, _13_ , e1631
  * Bao et al. 2018 Bao, J. L.; Gagliardi, L.; Truhlar, D. G. Self-interaction error in density functional theory: An appraisal. _The journal of physical chemistry letters_ 2018, _9_ , 2353–2358
  * Chołuj et al. 2018 Chołuj, M.; Kozłowska, J.; Bartkowiak, W. Benchmarking DFT methods on linear and nonlinear electric properties of spatially confined molecules. _International Journal of Quantum Chemistry_ 2018, _118_ , e25666
  * Stoll and Savin 1985 Stoll, H.; Savin, A. In _Density functional methods in physics_ ; da Providencia, R. M. D. . J., Ed.; Plenum, New York, 1985
  * Leininger et al. 1997 Leininger, T.; Stoll, H.; Werner, H.-J.; Savin, A. Combining long-range configuration interaction with short-range density functionals. _Chemical physics letters_ 1997, _275_ , 151–160
  * Savin 1996 Savin, A. In _In Recent Developments of Modern Density Functional Theory; Seminario, JM, Ed_ ; Seminario, J. M., Ed.; Elsevier: Amsterdam, 1996
  * Iikura et al. 2001 Iikura, H.; Tsuneda, T.; Yanai, T.; Hirao, K. A long-range correction scheme for generalized-gradient-approximation exchange functionals. _The Journal of Chemical Physics_ 2001, _115_ , 3540–3544
  * Yanai et al. 2004 Yanai, T.; Tew, D. P.; Handy, N. C. A new hybrid exchange–correlation functional using the Coulomb-attenuating method (CAM-B3LYP). _Chem. Phys. Lett._ 2004, _393_ , 51 – 57
  * 31 Built using the most popular OT technique that consists of enforcing the analogous of Koopmans’ theorem for DFT koopmans1934zuordnung, janak_prb1978.
  * Goerigk and Grimme 2014 Goerigk, L.; Grimme, S. Double-hybrid density functionals. _Wiley Interdisciplinary Reviews: Computational Molecular Science_ 2014, _4_ , 576–600
  * Bremond et al. 2016 Bremond, E.; Ciofini, I.; Sancho-García, J. C.; Adamo, C. Nonempirical double-hybrid functionals: An effective tool for chemists. _Accounts of chemical research_ 2016, _49_ , 1503–1513
  * Grimme 2006 Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. _J. Comput. Chem._ 2006, _27_ , 1787–1799
  * Brémond and Adamo 2011 Brémond, E.; Adamo, C. Seeking for parameter-free double-hybrid functionals: The PBE0-DH model. _The Journal of chemical physics_ 2011, _135_ , 024106
  * Brémond et al. 2014 Brémond, É.; Sancho-García, J. C.; Pérez-Jiménez, Á. J.; Adamo, C. Communication: Double-hybrid functionals from adiabatic-connection: The QIDH model. _The Journal of Chemical Physics_ 2014, _141_ , 031101
  * Langreth and Perdew 1977 Langreth, D. C.; Perdew, J. P. Exchange-correlation energy of a metallic surface: Wave-vector analysis. _Phys. Rev. B_ 1977, _15_ , 2884–2901
  * Furche and Van Voorhis 2005 Furche, F.; Van Voorhis, T. Fluctuation-dissipation theorem density-functional theory. _The Journal of chemical physics_ 2005, _122_ , 164106
  * Fuchs and Gonze 2002 Fuchs, M.; Gonze, X. Accurate density functionals: Approaches using the adiabatic-connection fluctuation-dissipation theorem. _Phys. Rev. B_ 2002, _65_ , 235109
  * Dobson and Wang 1999 Dobson, J. F.; Wang, J. Successful Test of a Seamless van der Waals Density Functional. _Phys. Rev. Lett._ 1999, _82_ , 2123–2126
  * Harl and Kresse 2008 Harl, J.; Kresse, G. Cohesive energy curves for noble gas solids calculated by adiabatic connection fluctuation-dissipation theory. _Phys. Rev. B_ 2008, _77_ , 045136
  * Schimka et al. 2010 Schimka, L.; Harl, J.; Stroppa, A.; Grueneis, A.; Marsman, M.; Mittendorfer, F.; Kresse, G. Accurate surface and adsorption energies from many-body perturbation theory. _Nat. Mater._ 2010, _9_ , 741–744
  * Lebègue et al. 2010 Lebègue, S.; Harl, J.; Gould, T.; Ángyán, J. G.; Kresse, G.; Dobson, J. F. Cohesive Properties and Asymptotics of the Dispersion Interaction in Graphite by the Random Phase Approximation. _Phys. Rev. Lett._ 2010, _105_ , 196401
  * Ren et al. 2012 Ren, X.; Rinke, P.; Blum, V.; Wieferink, J.; Tkatchenko, A.; Sanfilippo, A.; Reuter, K.; Scheffler, M. Resolution-of-identity approach to Hartree–Fock, hybrid density functionals, RPA, MP2 and GW with numeric atom-centered orbital basis functions. _New Journal of Physics_ 2012, _14_ , 053020
  * Ren et al. 2009 Ren, X.; Rinke, P.; Scheffler, M. Exploring the random phase approximation: Application to CO adsorbed on Cu(111). _Phys. Rev. B_ 2009, _80_ , 045402
  * Bruneval 2012 Bruneval, F. Range-Separated Approach to the RPA Correlation Applied to the van der Waals Bond and to Diffusion of Defects. _Phys. Rev. Lett._ 2012, _108_ , 256403
  * Bruneval et al. 2021 Bruneval, F.; Rodriguez-Mayorga, M.; Rinke, P.; Dvorak, M. Improved One-Shot Total Energies from the Linearized GW Density Matrix. _Journal of Chemical Theory and Computation_ 2021, _17_ , 2126–2136, PMID: 33705127
  * Rodríguez-Mayorga et al. 2021 Rodríguez-Mayorga, M.; Mitxelena, I.; Bruneval, F.; Piris, M. Coupling Natural Orbital Functional Theory and Many-Body Perturbation Theory by Using Nondynamically Correlated Canonical Orbitals. _Journal of chemical theory and computation_ 2021, _17_ , 7562–7574
  * Cohen and Handy 2001 Cohen, A. J.; Handy, N. C. Dynamic correlation. _Molecular Physics_ 2001, _99_ , 607–615
  * Rojas et al. 1995 Rojas, H. N.; Godby, R. W.; Needs, R. J. Space-Time Method for Ab Initio Calculations of Self-Energies and Dielectric Response Functions of Solids. _Phys. Rev. Lett._ 1995, _74_ , 1827–1830
  * Foerster et al. 2011 Foerster, D.; Koval, P.; Sánchez-Portal, D. An O(N3) implementation of Hedinś GW approximation for molecules. _J. Chem. Phys._ 2011, _135_ , 074105
  * Kaltak and Kresse 2020 Kaltak, M.; Kresse, G. Minimax isometry method: A compressive sensing approach for Matsubara summation in many-body perturbation theory. _Phys. Rev. B_ 2020, _101_ , 205145
  * Wilhelm et al. 2018 Wilhelm, J.; Golze, D.; Talirz, L.; Hutter, J.; Pignedoli, C. A. Toward GW Calculations on Thousands of Atoms. _JPCL_ 2018, _9_ , 306–312
  * 54 That may include a certain amount of exact exchange.
  * Bruneval et al. 2016 Bruneval, F.; Rangel, T.; Hamed, S. M.; Shao, M.; Yang, C.; Neaton, J. B. Molgw 1: Many-body perturbation theory software for atoms, molecules, and clusters. _Comput. Phys. Commun._ 2016, _208_ , 149 – 161
  * 56 Bruneval, F.; et al., https://github.com/bruneval/molgw
  * Rutishauser 1963 Rutishauser, H. Ausdehnung des Rombergschen Prinzips. _Numerische Mathematik_ 1963, _5_ , 48–54
  * Medveď et al. 2007 Medveď, M.; Stachová, M.; Jacquemin, D.; André, J.-M.; Perpete, E. A. A generalized Romberg differentiation procedure for calculation of hyperpolarizabilities. _Journal of Molecular Structure: THEOCHEM_ 2007, _847_ , 39–46
  * 59 The latest release of MOLGW code (version 3.2) incorporates these new features (download from https://github.com/bruneval/molgw/releases/tag/v3.2).
  * Dunning 1989 Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. _J. Chem. Phys._ 1989, _90_ , 1007–1023
  * Woon and Dunning 1993 Woon, D. E.; Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. III. The atoms aluminum through argon. _J. Chem. Phys._ 1993, _98_ , 1358–1371
  * Kendall et al. 1992 Kendall, R. A.; Dunning, T. H.; Harrison, R. J. Electron affinities of the first‐row atoms revisited. Systematic basis sets and wave functions. _J. Chem. Phys._ 1992, _96_ , 6796–6806
  * Hättig et al. 2006 Hättig, C.; Hellweg, A.; Köhn, A. Distributed memory parallel implementation of energies and gradients for second-order Møller–Plesset perturbation theory with the resolution-of-the-identity approximation. _Phys. Chem. Chem. Phys._ 2006, _8_ , 1159–1169
  * Weigend et al. 2002 Weigend, F.; Köhn, A.; Hättig, C. Efficient use of the correlation consistent basis sets in resolution of the identity MP2 calculations. _J. Chem. Phys._ 2002, _116_ , 3175–3183
  * 65 Molecular properties virtual lab, https://molprolab.com/data/
  * Bruneval 2019 Bruneval, F. Assessment of the Linearized GW Density Matrix for Molecules. _J. Chem. Theory Comput._ 2019, _15_ , 4069–4078
  * Denawi et al. 2023 Denawi, A. H.; Bruneval, F.; Torrent, M.; Rodríguez-Mayorga, M. GW density matrix in solids as an estimate of self-consistent GW total energies. _arXiv preprint arXiv:2304.13296_ 2023,
  * Perdew et al. 1996 Perdew, J. P.; Ernzerhof, M.; Burke, K. Rationale for mixing exact exchange with density functional approximations. _The Journal of chemical physics_ 1996, _105_ , 9982–9985
  * 69 Notice that for RPA@PBEh(0.0), the MAX% for $\alpha_{zz}$ is 38.47%, the MEAN% of $\gamma_{zzzz}$ is 96.85%, and the MAX% of $\gamma_{zzzz}$ is 898.91%; hence, these values lie out of the scale used.
  * Ren et al. 2013 Ren, X.; Rinke, P.; Scuseria, G. E.; Scheffler, M. Renormalized second-order perturbation theory for the electron correlation energy: Concept, implementation, and benchmarks. _Phys. Rev. B_ 2013, _88_ , 035120


# Refining perovskite structures to pair distribution function data using
collective Glazer modes as a basis

Sandra Helen Skjærvø  Department of Applied Physics and Applied Mathematics,
Columbia University, New York, NY, 10027, USA  Martin A. Karlsen  Department
of Physics, Chemistry and Pharmacy, University of Southern Denmark, 5230,
Odense M, Denmark  Riccardo Comin  Physics Department, Massachusetts Institute
of Technology, Cambridge, MA, 02139, USA  Simon J. L. Billinge  Department of
Applied Physics and Applied Mathematics, Columbia University, New York, NY,
10027, USA  Condensed Matter Physics and Materials Science Department,
Brookhaven National Laboratory, Upton, NY, 11973, USA

###### Abstract

Structural modelling of octahedral tilts in perovskites is typically done
using the symmetry constraints of the resulting space group. In most cases,
this introduces more degrees of freedom than those strictly necessary to
describe only the octahedral tilts. It can therefore be a challenge to
disentangle the octahedral tilts from other structural distortions such as
cation displacements and octahedral distortions. This paper reports on the
development of constraints for modelling pure octahedral tilts and implemented
the constraints in diffpy-CMI, a powerful package to analyse pair distribution
function (PDF) data. The program allows features in the PDF that come from
rigid tilts to be separated from non-rigid relaxations, provides an intuitive
picture of the tilting, and as it has many fewer refinable variables than the
unconstrained space-group fits, provides robust and stable refinements of the
tilt components. It further demonstrates the use of the model on the canonical
tilted perovskite CaTiO3 which has a known Glazer tilt system
$\alpha^{+}\beta^{-}\beta^{-}$. The Glazer model fits comparably to the
corresponding space group model $Pnma$ below $r=14$ Å and becomes
progressively worse than the space group model at higher $r$ due to non-rigid
distortions in the real material.

##  1 Introduction

Structural distortions in materials, such as those occurring during displacive
structural phase transitions, often involve collective displacements of groups
of atoms [Dove, 1997]. For example, in the perovskites, a material class with
nominal stoichiometry ABX3 (Fig. 1), collective distortions are known to cause
a host of structural phase transitions that lower the symmetry of the cubic
parent structure [Müller et al., 1968, Salje, 1990, Goodenough, 1955, Kwei et
al., 1993, Kwei et al., 1995].

Figure 1: Illustration of in-phase and out-of-phase tilt systems as viewed
down the tilt axis. The tilt systems shown here are
$\alpha^{0}\alpha^{0}\gamma^{+}$ (top) and $\alpha^{0}\alpha^{0}\gamma^{-}$
(bottom).

Distortions away from from the cubic archetype can involve deformations of the
octahedra, displacements of the B cations inside the octahedra and tilting of
the octahedra. The first two are typically caused by electronic instabilities,
while the latter is due to the relative sizes of the cations. For perovskites
with smaller A cations, the octahedra tilt to compress the structure around
them, essentially improving the bonding for the A cation. This geometric
effect is conveniently captured by the Goldschmidt tolerance factor
[Goldschmidt, 1926],

$t=\frac{r_{\mathrm{A}}+r_{\mathrm{X}}}{\sqrt{2}(r_{\mathrm{B}}+r_{\mathrm{X}})}$
(1)

where $r$ is an ionic radius and subscripts A, B and X denote the ion type.
For $t=1$, the perovskite crystallizes in the high-symmetry cubic structure,
while octahedral tilting is expected for a $t<1$ as it signifies that the A
site cation is too small to fill the void between the octahedra. In this paper
we will concentrate on the latter type of distortion.

Due to their corner-sharing geometry, the octahedra can tilt collectively in
several different patterns. By building macroscopic models of corner-shared
rigid octahedra, Glazer was able to describe all the 22 different patterns in
which the rigid octahedra could collectively tilt, and the resulting symmetry
space groups [Glazer, 1972, Glazer, 1975]. Later studies uncovered details
about these Glazer systems through group theory and geometric considerations
[Aleksandrov, 1976, O’Keeffe and Hyde, 1977, Woodward, 1997a, Woodward, 1997b,
Howard and Stokes, 1998].

Depending on the Glazer tilt pattern of a perovskite, the structure will have
different symmetry space group [Aleksandrov, 1976, O’Keeffe and Hyde, 1977,
Woodward, 1997a, Woodward, 1997b, Howard and Stokes, 1998]. Modelling the
structures of the these low-symmetry phases is therefore often done using the
symmetry-broken crystallographic models and constraining the allowed atomic
displacements to those imposed by the space-group symmetries. However, in
general, the symmetry space group allows for more displacive degrees of
freedom than those strictly needed to describe the tilting of the octahedra.
Using these models for fitting scattering data leads to structures where the
octahedra are distorted in a way that cannot be represented in terms of the
pure Glazer tilt patterns with rigid units even in the cases where the
octahedra are not geometrically required to distort [Howard and Stokes, 1998].

Here we explore a more direct approach to modeling collective rotations by
using algebraic expressions that link displacements of atoms in the Glazer
tilt systems. Going beyond purely symmetry constraints is surprisingly
challenging. Approximate Monte Carlo approaches have been attempted [Sartbaeva
et al., 2006, Sartbaeva et al., 2007], where atoms are tethered to rigid-unit
templates which do not distort, but are allowed to relax away from the
vertices. It has also been shown [Campbell et al., 2018] that for small
rotations a set of linear equations on top of symmetry mode analysis [Perez-
Mato et al., 2010] can identify collective modes in a system of connected
rigid units that do not (or hardly) distort the units. However, there is
currently no straightforward way of incorporating this information into a
refinement program for quantitative modeling of data in terms of this
collective mode basis.

Our approach of explicitly building the geometric constraint equations without
assumed symmetries has an advantage that it can be easily plugged into local
structure modeling schemes such as that used in the diffpy-CMI [Juhás et al.,
2015] program. The program works in the P1 space-group by design, allowing one
to introduce structural distortions by moving atoms at will. The approach
greatly reduces the number of refinable parameters in a physically meaningful
way and can help to build intuition about the structure and how it is likely
to distort. It also allows the user to test directly hypotheses about the
rigidity of the units or the type of tilting present in a sample without the
conceptual complexity of having to surf between space groups. This can give
new insight that might be lost otherwise. Here we demonstrate the use of our
code on the compound CaTiO3, the archetypal perovskite with a well-known
Glazer tilt pattern $\alpha^{+}\beta^{-}\beta^{-}$.

##  2 Glazer tilt definitions

The Glazer tilt systems, as laid out by Glazer in 1972 [Glazer, 1972],
describe the complete set of collective rotations allowed in a network of
corner-shared octahedra as found in perovskites111The tilt patterns described
by Glazer can all be described using a 2x2x2 (or smaller) supercell of the
cubic perovskite unit cell, and collective distortions requiring larger
supercells are unlikely., shown in Table 1.

Table 1: For each of the different Glazer tilt patterns we give the index as assigned by Glazer (Glazer, A. M. 1972; Glazer, A. M. 1975), the tilts given with Glazer notation and the space group symmetry of the resulting phase. Note that we have only included here the tilt systems that are symmetry nonequivalent. Tilt |  |
---|---|---
system |  Tilts |  Space group
23 |  $\alpha^{0}\alpha^{0}\alpha^{0}$ |  $Pm\bar{3}m$ (#221)
22 |  $\alpha^{0}\alpha^{0}\gamma^{-}$ |  $I4/mcm$ (#140)
21 |  $\alpha^{0}\alpha^{0}\gamma^{+}$ |  $P4/mbm$ (#127)
20 |  $\alpha^{0}\beta^{-}\beta^{-}$ |  $Imma$ (#74)
19 |  $\alpha^{0}\beta^{-}\gamma^{-}$ |  $C2/m$ (#12)
17 |  $\alpha^{0}\beta^{+}\gamma^{-}$ |  $Cmcm$ (#63)
16 |  $\alpha^{0}\beta^{+}\beta^{+}$ |  $I4/mmm$ (#139)
14 |  $\alpha^{-}\alpha^{-}\alpha^{-}$ |  $R\bar{3}c$ (#167)
13 |  $\alpha^{-}\beta^{-}\beta^{-}$ |  $C2/c$ (#15)
12 |  $\alpha^{-}\beta^{-}\gamma^{-}$ |  $P\bar{1}$ (#2)
10 |  $\alpha^{+}\beta^{-}\beta^{-}$ |  $Pnma$ (#62)
8 |  $\alpha^{+}\beta^{-}\gamma^{-}$ |  $P2_{1}/m$ (#11)
5 |  $\alpha^{+}\alpha^{+}\gamma^{-}$ |  $P4_{2}/nmc$ (#137)
3 |  $\alpha^{+}\alpha^{+}\alpha^{+}$ |  $Im\bar{3}$ (#204)
1 |  $\alpha^{+}\beta^{+}\gamma^{+}$ |  $Immm$ (#71)

For the sake of clarity, we introduce the naming scheme here. An octahedron
can be tilted around one, two or all three of the cartesian axes, $x$, $y$ and
$z$. The nature of each rotation is indicated by three Greek letters with
superscripts, where the first letter denotes the rotation around $x$, the
second around $y$ and the third around $z$. Repeating letters (e.g.
$\alpha^{+}\alpha^{+}\alpha^{+}$) indicate that the amplitudes around the
specific axes are the same, while different letters (e.g.
$\alpha^{+}\beta^{+}\gamma^{+}$) indicate that the tilts differ in amplitude
around the different axes.

The superscripts can take the value 0, + or - to indicate a zero-tilt
amplitude or a non-zero amplitude with tilts in adjacent layers along the tilt
axis being either in-phase (+) or out-of-phase (-). For example, the tilt
pattern $\alpha^{0}\alpha^{0}\gamma^{+}$ has no tilt around the $x$ and $y$
axes and a non-zero tilt around the $z$ axis. Because of the connectivity of
the octahedra at their corners, neighbouring octahedra in the plane
perpendicular to the tilt axis rotate in the opposite direction to the central
octahedron, leading to a doubling of the unit cell in that plane. In the
example of $\alpha^{0}\alpha^{0}\gamma^{+}$, the unit cell is therefore
doubled in the $ab$ plane, but not along the $z$ axis. On the other hand, an
out-of-phase tilt, for example along the $z$ axis in the pattern
$\alpha^{0}\alpha^{0}\gamma^{-}$, will double the unit cell also along the
tilt axis. Fig. 1 illustrates the difference between the in-phase and out-of-
phase tilt pattern of the $\alpha^{0}\alpha^{0}\gamma^{+}$ and
$\alpha^{0}\alpha^{0}\gamma^{-}$ tilt systems, as viewed down the tilt axis.

##  3 Approach

Here we describe the method for building constrained Glazer tilt pattern
models. The code may be found at https://github.com/sandraskj/glazer_fitting.

Models are built using the diffpy-CMI program [Juhás et al., 2015], which has
powerful and flexible methods for specifying constraints between model
parameters. This allows, in principle, large numbers of parameters to be
expressed in terms of a much smaller number of variables from analytic or
numerical expressions. We first generate the constraints as symbolic
expressions relating multiple atoms’ fractional coordinates. These expressions
are then captured into the diffpy-CMI constraint handling interface.

For all the Glazer tilt systems listed in Table 1 the constraints have been
constructed such that the shortest B-X distances are all kept rigid. Since,
for most of the systems, there is a small coupling between the rotation modes
around the three axes, these constraints will lead to a small octahedral
distortion and octahedral angles deviating slightly from 90∘, so the tilts are
not strictly rigid. However, the tilt equations result in almost rigid
octahedral tilting, where the collective modes may be described just in terms
of tilt angles around each axis which are the only refinable parameters for
the modes when fitting to data, in addition to the cubic lattice parameter.

The collective octahedral rotations do not include A-site ion structural
parameters. Although the A atoms are not directly part of the octahedral
tilting network, their positions are still refined, as they do respond to the
tilts by displacing. We chose to constrain the A cation displacements in such
a way as to respect the expected symmetry of the tilted structure, which in
the case of CaTiO3 is the space group $Pnma$.

Activating tilt modes leads to well-defined reductions in the lattice
parameters, and therefore for a full description, we need to find the
appropriate scaling parameters expressed in terms of the Glazer tilt
amplitudes and the lattice parameter of the cubic parent structure. We start
with the interatomic vectors from the B atom at the origin to its three unique
X neighbors in the octahedron, $\mathbf{r}_{\mathrm{B-X1}}$,
$\mathbf{r}_{\mathrm{B-X2}}$ and $\mathbf{r}_{\mathrm{B-X3}}$,

$\displaystyle|\mathbf{r}_{\mathrm{B-X1}}|$
$\displaystyle=\sqrt{\left(a^{\prime}\frac{1}{4}\right)^{2}+\left(b^{\prime}\cdot
y_{\mathrm{X1}}\right)^{2}+\left(c^{\prime}\cdot
z_{\mathrm{X1}}\right)^{2}}\,,$ (2)
$\displaystyle|\mathbf{r}_{\mathrm{B-X2}}|$
$\displaystyle=\sqrt{\left(a^{\prime}\cdot
x_{\mathrm{X2}}\right)^{2}+\left(b^{\prime}\frac{1}{4}\right)^{2}+\left(c^{\prime}\cdot
z_{\mathrm{X2}}\right)^{2}}\,,$ (3)
$\displaystyle|\mathbf{r}_{\mathrm{B-X3}}|$
$\displaystyle=\sqrt{\left(a^{\prime}\cdot
x_{\mathrm{X3}}\right)^{2}+\left(b^{\prime}\cdot
y_{\mathrm{X3}}\right)^{2}+\left(c^{\prime}\frac{1}{4}\right)^{2}}\,,$ (4)

where $a^{\prime}$, $b^{\prime}$ and $c^{\prime}$ are the lattice parameters
of the distorted supercell for a given set of Glazer tilts. Keep in mind that
the fractional coordinates $y_{\mathrm{X1}}$, $z_{\mathrm{X2}}$ etc., are all
expressions containing the Glazer tilt variables and the lattice parameter of
the cubic parent structure $a_{h}$. Next, we set each of the bond lengths to
be a quarter of the parent unit cell $a_{h}$, i.e., for X1,

$\displaystyle|\mathbf{r}_{\mathrm{B-X1}_{x}}|=|\mathbf{r}_{\mathrm{B-X1}_{y}}|=|\mathbf{r}_{\mathrm{B-X1}_{z}}|=a_{h}/4\,,$
(5)

where $\mathbf{r}_{\textrm{B-X1}_{z}}$ is the $z$ component of
$\mathbf{r}_{\textrm{B-X1}}$. Since the rotations are assumed to be rotations
of rigid octahedra, these lengths will not change after the rotation. This
allows us to relate the lattice parameter of the Glazer tilt distorted
supercell to the ones of the cubic parent cell through scaling parameters
$s_{a}$, $s_{b}$ and $s_{c}$,

$\displaystyle
a^{\prime}=\frac{2a_{h}}{s_{a}}\,,\;b^{\prime}=\frac{2a_{h}}{s_{b}}\,,\;c^{\prime}=\frac{2a_{h}}{s_{c}}.$
(6)

Substituting for $a^{\prime}$, $b^{\prime}$ and $c^{\prime}$ in Eqs. 2-4, we
get a set of three equations,

$\displaystyle|\mathbf{r}_{\mathrm{B-X1}}|=\sqrt{\frac{1}{16s_{a}}+\left(\frac{y_{\mathrm{X1}}}{s_{b}}\right)^{2}+\left(\frac{z_{\mathrm{X1}}}{s_{c}}\right)^{2}}=a_{h}/4\,,$
(7)
$\displaystyle|\mathbf{r}_{\mathrm{B-X2}}|=\sqrt{\left(\frac{x_{\mathrm{X2}}}{s_{a}}\right)^{2}+\frac{1}{16s_{b}}+\left(\frac{z_{\mathrm{X2}}}{s_{c}}\right)^{2}}=a_{h}/4\,,$
(8)
$\displaystyle|\mathbf{r}_{\mathrm{B-X3}}|=\sqrt{\left(\frac{x_{\mathrm{X3}}}{s_{a}}\right)^{2}+\left(\frac{y_{\mathrm{X3}}}{s_{b}}\right)^{2}+\frac{1}{16s_{c}}}=a_{h}/4\,,$
(9)

that can be solved for $s_{a}$, $s_{b}$ and $s_{c}$,

$\displaystyle s_{a}$
$\displaystyle=\sqrt{\frac{1-256x_{\mathrm{X2}}^{2}\left(y_{\mathrm{X1}}^{2}-16y_{\mathrm{X3}}^{2}z_{\mathrm{X1}}^{2}\right)-256\left(y_{\mathrm{X3}}^{2}z_{\mathrm{X2}}^{2}+x_{\mathrm{X3}}^{2}\left(z_{\mathrm{X1}}^{2}-16y_{\mathrm{X1}}^{2}z_{\mathrm{X2}}^{2}\right)\right)}{1-16y_{\mathrm{X1}}^{2}+16\left(-1+16y_{\mathrm{X3}}^{2}\right)z_{\mathrm{X1}}^{2}+256(y_{\mathrm{X1}}^{2}-y_{\mathrm{X3}}^{2})z_{\mathrm{X2}}^{2}}}\,,$
(10) $\displaystyle s_{b}$
$\displaystyle=\sqrt{\frac{1-256x_{\mathrm{X2}}^{2}\left(y_{\mathrm{X1}}^{2}-16y_{\mathrm{X3}}^{2}z_{\mathrm{X1}}^{2}\right)-256\left(y_{\mathrm{X3}}^{2}z_{\mathrm{X2}}^{2}+x_{\mathrm{X3}}^{2}\left(z_{\mathrm{X1}}^{2}-16y_{\mathrm{X1}}^{2}z_{\mathrm{X2}}^{2}\right)\right)}{1+16x_{\mathrm{X2}}^{2}\left(-1+16z_{\mathrm{X1}}^{2}\right)-16z_{\mathrm{X2}}^{2}+256x_{\mathrm{X3}}^{2}\left(-z_{\mathrm{X1}}^{2}+z_{\mathrm{X2}}^{2}\right)}}\,,$
(11) $\displaystyle s_{c}$
$\displaystyle=\sqrt{\frac{-1-256x_{\mathrm{X2}}^{2}\left(y_{\mathrm{X1}}^{2}-16y_{\mathrm{X3}}^{2}z_{\mathrm{X1}}^{2}\right)+256\left(y_{\mathrm{X3}}^{2}z_{\mathrm{X2}}^{2}+x_{\mathrm{X3}}^{2}\left(z_{\mathrm{X1}}^{2}+16y_{\mathrm{X1}}^{2}z_{\mathrm{X2}}^{2}\right)\right)}{-1+x_{\mathrm{X3}}^{2}(16-256y_{\mathrm{X1}}^{2})+256x_{\mathrm{X2}}^{2}(y_{\mathrm{X1}}^{2}-y_{\mathrm{X3}}^{2})+16y_{\mathrm{X3}}^{2}}}\,.$
(12)

Setting these as constraints in the refinement allows the unit cell to change
size according to the tilt amplitude without introducing any extra refinable
parameters.

We present the full constraints for Glazer system 10, as generated from the
code, in the project code repository on GitHub
(https://github.com/sandraskj/glazer_fitting). In the GitHub repository we
also provide the code that generates the constraints for all the Glazer
systems and brief instructions for how the reader can download them and how to
set it up for their own refinements using diffpy-CMI.

##  4 Experimental measurements

To obtain experimental pair distribution functions for CaTiO3 measurements
were carried out at the 28-ID-2 (XPD) beamline at the NSLS-II at Brookhaven
National Laboratory on a commercially purchased powder sample of CaTiO3 (Strem
Chemicals Inc, CAS 12049-50-2). A 2D Perkin Elmer amorphous silicon detector
was placed 380 mm behind the sample, which was loaded in a 0.5 mm glass
capillary. The wavelength of the incident x-rays was 0.16635 Å. Data were
collected at 200 K for 60 s in a flowing nitrogen cryostream.

The data were processed using standard methods [Egami and Billinge, 2012]. The
instrument geometry was calibrated using data from a fine powdered Ni powder
using pyFAI [Kieffer et al., 2020]. 2D diffraction patterns were processed by
applying masks to remove the beam stop as well as outlier saturated and dead
pixels using a home-written automasking protocol. After correction for
polarization effects they were integrated azimuthally along circles of
constant scattering vector magnitude, $Q$, also using pyFAI. The background
signal from an empty glass capillary was subtracted and the data were
normalized and corrected to obtain the reduced total scattering structure
function, $F(Q)$, which was Fourier transformed to obtain the PDF. This was
done using PDFgetX3 [Juhás et al., 2013]. The maximum range of data used in
the Fourier transform was $Q_{max}=23.6$ Å-1

##  5 Results

Our initial tests of the approach are carried out on simulated data from the
known ground-state structure of CaTiO3.

Figure 2: Comparison of the fit of the 14 three-tilt Glazer systems to a
calculated PDF of CaTiO3 with octahedral rotations but without Ca
displacements.

The structure was created in Glazer system no. 10 with an in-phase tilt around
one axis of $\alpha=9^{\circ}$ and out-of-phase tilt around the other two axes
of $\beta=10^{\circ}$. For simplicity the simulated structure had no
displacement of the Ca atoms away from their cubic positions. The isotropic
atomic displacement parameters (ADPs) for all the ions were set to
$U_{iso}(\mathrm{Ca})=0.0030$ Å2, $U_{iso}(\mathrm{Ti})=0.0046$ Å2, and
$U_{iso}(\mathrm{O})=0.011$ Å2, similar to those obtained from fitting CaTiO3
with a conventional $Pnma$ model.

The PDF was calculated from the structure using diffpy-CMI [Juhás et al.,
2015], with damping and broadening parameters set to values 0.029 Å-1 and
0.010 Å-1, respectively, obtained from the calibration sample in our
measurement and $Q_{max}=23.6$ Å-1, the same value as we used for the
experimental data.

We then fit constrained Glazer models from each of the 22 Glazer tilt patterns
to the data from $1.6<r<15$ Å to see how well the constrained refinements were
working. The starting values for the tilt amplitudes in the refinements models
were set to values that were roughly 70% of the true values in the structure
for the calculated data set.

The fits with the one-tilt and two-tilt models were poor in most cases, while
all the three-tilt systems gave fit residuals below 10%. Tilt system 10 (the
ground-truth result) is one of the three-tilt systems so this gives confidence
that the approach can easily differentiate the presence or absence of tilts.
However, within the subset of three-tilt systems, different families of tilt
combinations can be found which refine to significantly different $R_{w}$
values, as shown in Fig. 2. Interestingly, the fits can differentiate cases
that have $+++$, $++-$, $+--$ and $---$ tilts, but within those families it
cannot distinguish between different Glazer systems. This may be because the
tilt amplitudes we chose for the test, coming from the observed values in
CaTiO3, are close to each other.

The best overall fit was found for Glazer system 10, the correct one, as well
as Glazer system 8 that has the same tilt pattern, but with an extra degree of
freedom that allows the out-of-phase tilts to be of different amplitudes. This
shows that the collective mode refinements are working in diffpy-CMI.

Table 2: Comparison of parameters from the space group and Glazer model refinements over the $r$ range 1.6-50 Å. The parameters typeset in bold were variables that were explicitly refined, while the parameters typeset in italic were extracted or calculated post-refinement from the refined values. Two values are given each for $\alpha$ and $\beta$ in the space group model because different octahedra tilt by different amounts. $Pnma$ space group model | Glazer model
---|---
variable | value | variable | value
scale | 0.18 | scale | 0.17
delta1 | 1.03 | delta1 | 2.47
|  | $\mathbf{a_{h}}$ | 3.907
$\mathbf{a}$ | 5.428 | $a$ | 5.402
$\mathbf{b}$ | 7.620 | $b$ | 7.594
$\mathbf{c}$ | 5.366 | $c$ | 5.402
$\mathbf{x_{Ca}}$ | 0.0357 | $\mathbf{x_{Ca}}$ | 0.0216
$\mathbf{z_{Ca}}$ | 0.0031 | $\mathbf{y_{Ca}}$ | 0.0069
$\alpha$ | 9.6∘, 10.1∘ | $\boldsymbol{\alpha}$ | 7.6∘
$\beta$ | 7.1∘, 10.6∘ | $\boldsymbol{\beta}$ | 9.7∘
$\mathbf{x_{O1}}$ | 0.2059 |  |
$\mathbf{y_{O1}}$ | 0.0335 |  |
$\mathbf{z_{O1}}$ | 0.2073 |  |
$\mathbf{x_{O2}}$ | 0.0155 |  |
$\mathbf{z_{O2}}$ | 0.5784 |  |
$\mathbf{U_{iso}}$(Ca) | 0.005 |  $\mathbf{U_{iso}}$(Ca) | 0.004
$\mathbf{U_{iso}}$(Ti) | 0.003 |  $\mathbf{U_{iso}}$(Ti) | 0.004
$\mathbf{U_{iso}}$(O) | 0.010 |  $\mathbf{U_{iso}}$(O) | 0.011
$R_{w}$ | 0.087 | $R_{w}$ | 0.245

Next, we performed refinements on an experimental data set of CaTiO3. We
performed the refinements with two models: One using our formulation based on
Glazer tilt system 10 and for comparison, a model with constraints consistent
with the crystallographic space-group $Pnma$, which does not impose rigid
tilts. The Ca sites were constrained the same way in both models, according to
the space group symmetry of $Pnma$. The space group model has 10 structural
degrees of freedom, while the Glazer model has only 5. The variables,
including explicitly refined as well as post-calculated ones, and their values
after refinement over 1.6-50 Å are listed in Table 2.

Figure 3: Plots of measured (blue) and best-fit (red) PDFs of CaTiO3 with the
difference curve plotted in green offset below over the $r$ range 1.6-50Å. The
model for the best-fit PDF is from (a) the constrained Glazer tilt model in
Glazer system 10 and (b) allowing all the structural degrees of freedom from
the $Pnma$ structural model.

Fitting both models over this wide $r$ range (Fig. 3), we can see that the
space group model gives a significantly smaller fit residual (space group
gives $R_{w}$ = 0.087 while the Glazer model gives $R_{w}$ = 0.245), which is
not surprising given its larger number of refinable variables. Comparing the
refined structural parameters from the two models, we see that all, except for
the tilt angles and the lattice parameters are in quite good agreement as
shown by comparing the values in Table 2. The information of interest to us is
the presence and amplitude of rigid Glazer tilt modes. For the Glazer model
these are a direct output of the program. For comparison, we also calculate
the tilt angles from the space group by finding the angle that each vector
between opposite pairs of oxygen atoms on an octahedron makes with the
pseuodocubic axes. Previous studies have calculated the angles using the
oxygen positions similarly to what we do here, but in such a way that gives
the average over two octahedra in order to to end up with multiple values for
each Glazer component [Kennedy et al., 1999, Yashima and Ali, 2009]. This is
necessary if you want a single value for each tilt component as the space
group model does not keep the octahedral bond lengths rigid, and the tilts
therefore vary depending on the bond chosen for evaluating it. For the sake of
comparing the performance of our Glazer model to the space group model, we
choose here to not present such average tilt angles, but rather present the
tilts are they present themselves on different octahedra. This will highlight
the difference in rigidity and robustness of the two models.

The Glazer model results in values of in-phase tilt $\alpha=7.6^{\circ}$ and
out-of-phase tilt $\beta=9.7^{\circ}$. The space group model gives the values
$\alpha=9.6^{\circ}$ and 10.1∘ and $\beta=7.1^{\circ}$ and 10.6∘. The average
value of the in-phase $\alpha$ tilt is higher in the space-group model by
almost 2∘ than in the Glazer model, and that of the out-of-phase $\beta$ tilt
is lower by about 1∘. In the case of the $\beta$ the two different octahedra
in the space group model are quite different and actually straddle the value
obtained in the Glazer fit.

Addressing the difference in fit quality, we believe that a significant
contribution to the poorer fit is due to the tighter constraints on the
lattice parameters in the Glazer fits. The space group model is orthorhombic,
with three different lattice parameters. For the Glazer model only one lattice
parameter variable is refined. The tilts then result in a change in shape of
the unit from cubic, as described above. We note that the particular tilt
pattern in this Glazer mode results in a tetragonal, not orthorhombic, unit
cell (Table 2). The $a=c$ parameters for the Glazer model lie in value between
those of the space-group model, but are not able to separate into short and
long values allowed by the orthorhombic crystallographic model due to the
Glazer model constraints, whereas clearly structural relaxations beyond the
rigid tilts are present in the actual material that prefer this.

Figure 4: Comparison of the a) fit residual $R_{w}$, the octahedral tilt
amplitudes b) $\alpha$ and c) $\beta$, d) the total Ca displacements
($\delta_{Ca}$) e) and the $U_{\mathrm{iso}}$ values from boxcar fits with the
space group model and the Glazer model of CaTiO3 at 200 K. The $r$ range (or
’the box’) was set to 8 Å and incrementally shifted to higher $r$ values in
steps of 2 Å. The labels on the $x$-axis correspond to the highest value in
the box, $r_{max}$. The dotted lines represent the values obtained from a fit
over the 1.6-50 Å range. We note that for the space group model, the tilt
angles $\alpha$ and $\beta$ differ depending on which octahedra were used to
calculate them, and such are represented by two different lines.

If the difference in $R_{w}$ between the two models in the wide-range fits is
due to the difference in model rigidity, the models would be expected to
perform more comparably when fitting only the most local structure, and for
the Glazer model to perform worse at higher values of $r$. The Glazer model
only allows for the degrees of freedom that are strictly necessary for the
tilt pattern $\alpha^{+}\beta^{-}\beta^{-}$, and comparing the two fits at
different length scales therefore allows us to separate contributions to the
PDF signal that come from rigid tilts and additional non-rigid relaxations. It
is also interesting to consider if the refined values of the Glazer tilts
varies with the $r$-range that is fit over, as might be the case if the tilts
become damped with increasing-$r$. We therefore performed a series of fits
where an $r$-range of a fixed size (referred to as a box) is shifted
incrementally up to higher values, an approach we call a ‘boxcar’ fit. The
$r$-dependence of the refined variables are shown in Fig. 4.

As evident in Fig. 4(b) and (c) the values of the tilt amplitudes vary more
smoothly in the Glazer model than in the space-group model indicating that
refinement of these variables is more stable in the more highly constrained
Glazer fits. Also, whilst the $\alpha$ tilt is fairly $r$-independent, there
is a marked tendency for the $\beta$ tilt to decrease with increasing-$r$ in
the Glazer fit. This would be expected if there is a loss in the coherence of
the tilts with increasing-$r$ due to a non-rigidity. This suggests that in
cases where local tilts survive in a material but are not present globally,
the range of coherence of the collective motions may be measured by this
approach. We also see a similar trend in the total displacement of Ca from
cubic positions ($\delta_{Ca}$, panel d), with the Glazer model trending
downwards, while the space group model stays at the same value throughout the
$r$ range.

Figure 5: Plots of measured (blue) and best-fit (red) PDFs of CaTiO3 with the
difference curve plotted in green offset below over an $r$ range of 1.6-14 Å.
The model for the best-fit PDF is from (a) the constrained Glazer tilt model
in Glazer system 10 and (b) allowing all the structural degrees of freedom
from the $Pnma$ structural model.

Since we believe that a non-rigidity in the tilts explains the difference in
behavior of the two fits, a comparison for the fits over a much narrower range
$r=1.6-14$ Å we might expect much better agreement between the Glazer and
space-group models. This is indeed found in the comparable $R_{w}$ values at
low-$r$ (Fig. 4(a)). We show the fits over this low-$r$ region on an expanded
scale in Fig. 5. This shows that the most local structure is sufficiently
rigid, with octahedral tilts being the predominant structural distortion, such
that it can be well represented by the Glazer model.

We note that for the case we studied, CaTiO3 at room temperature, the tilts
are long-range ordered and so are expected to persist over large distances,
asymptotically approaching the crystallographic values. This kind of boxcar
analysis can be expected to be more interesting in materials where no tilts
are observed in the average structure but are observed locally [Skjærvø et
al., 2019, Bozin et al., 2019, Koch et al., 2021, Yang et al., 2020, Wang et
al., 2020, Senn et al., 2016].

##  6 Conclusions

We have developed sets of constraint equations that explicitly model
octahedral tilts (Glazer tilts) in perovskites. The model allows refinements
of collective atomic motions by geometrically connecting atoms in the lattice
allowing rigid rotations to be modeled directly. We have implemented the
constraints directly in the PDF modeling program diffpy-CMI.

We have demostrated the use of our code on the canonical tilted perovskite
system CaTiO3, which has a known long-range ordered Glazer tilt system
$\alpha^{+}\beta^{-}\beta^{-}$. We found that our Glazer model fits comparably
to the known space group model $Pnma$ below $r=14$ Å. We further saw that the
Glazer model performed progressively worse at higher $r$, due to the rigidity
of the model. In this case the rigid tilts alone broke the cubic symmetry to
tetragonal, whereas the observed symmetry is orthorhombic, which explains the
discrepancy in the fit residuals. Presumably, non-rigid relaxations and
relaxations of atoms not involved in the tilts are responsible for the
additional reduction in symmetry.

The use of our Glazer model could be used to study a wide range of perovskite
systems to better understand whether their structure is well explained in
terms of pure octahedral rotations, how the rotations vary with parameters
such as temperature and pressure, and what additional structural relaxations
are needed to explain the structure beyond the simple picture of octahedral
rotations. The highly constrained fits can be expected to give stable
refinements even when data quality is limited, for example, from small
nanoparticles or powders in a diamond anvil cell. The work also highlights the
strengths and limitations of the geometric approach in building rigid body
constraints.

##  7 Acknowledgements

We acknowledge Daniel E. DeRosha and Jonathan Owen for help in preparing the
CaTiO3 sample, and Milinda Abeykoon and Gihan Kwon for help with collecting
data at the 28-ID-1 (PDF) beamline at the NSLS-II. This work was supported by
the U.S. National Science Foundation through grant DMREF-1922234. M.A.K.
acknowledges support from the Carlsberg Foundation (grant no. CF17-0823). Use
of the National Synchrotron Light Source II, Brookhaven National Laboratory,
was supported by the US Department of Energy, Office of Science, Office of
Basic Energy Sciences, under Contract No. DE-SC0012704.

## References

  * [Aleksandrov, 1976] Aleksandrov, K. S. (1976).  The sequences of structural phase transitions in perovskites.  Ferroelectrics, 14(1):801–805.
  * [Bozin et al., 2019] Bozin, E. S., Yin, W. G., Koch, R. J., Abeykoon, M., Hor, Y. S., Zheng, H., Lei, H. C., Petrovic, C., Mitchell, J. F., and Billinge, S. J. L. (2019).  Local orbital degeneracy lifting as a precursor to an orbital-selective Peierls transition.  Nature Communications, 10(1):1–7.
  * [Campbell et al., 2018] Campbell, B. J., Averett, T., Yost, C., Machlus, S., Stokes, H. T., Whittle, T. A., Schmid, S., Howard, C. J., and IUCr (2018).  An algebraic approach to cooperative rotations in networks of interconnected rigid units.  Acta Crystallographica Section A: Foundations and Advances, 74:a286–a286.
  * [Dove, 1997] Dove, M. T. (1997).  Theory of displacive phase transitions in minerals.  American Mineralogist, 82(3):213–244.
  * [Egami and Billinge, 2012] Egami, T. and Billinge, S. J. L. (2012).  Underneath the Bragg peaks: structural analysis of complex materials.  Number 16 in Pergamon materials series. Elsevier, Amsterdam, 2nd edition.
  * [Glazer, 1972] Glazer, A. M. (1972).  The classification of tilted octahedra in perovskites.  Acta Crystallographica Section B, 28(11):3384–3392.
  * [Glazer, 1975] Glazer, A. M. (1975).  Simple ways of determining perovskite structures.  Acta Crystallographica Section A: Crystal Physics, Diffraction, Theoretical and General Crystallography, 31(6):756–762.
  * [Goldschmidt, 1926] Goldschmidt, V. M. (1926).  Die Gesetze der Krystallochemie.  Naturwissenschaften, 14(21):477–485.
  * [Goodenough, 1955] Goodenough, J. B. (1955).  Theory of the Role of Covalence in the Perovskite-Type Manganites La, m(II)MnO3.  Physical Review, 100(2):564–573.
  * [Howard and Stokes, 1998] Howard, C. J. and Stokes, H. T. (1998).  Group-Theoretical Analysis of Octahedral Tilting in Perovskites.  Acta Crystallographica Section B, 54(6):782–789.
  * [Juhás et al., 2013] Juhás, P., Davis, T., Farrow, C. L., and Billinge, S. J. L. (2013).  PDFgetX3: a rapid and highly automatable program for processing powder diffraction data into total scattering pair distribution functions.  Journal of Applied Crystallography, 46(2):560–566.
  * [Juhás et al., 2015] Juhás, P., Farrow, C., Yang, X., Knox, K., and Billinge, S. J. L. (2015).  Complex modeling: a strategy and software program for combining multiple information sources to solve ill posed structure and nanostructure inverse problems.  Acta Crystallographica Section A, 71(6):562–568.
  * [Kennedy et al., 1999] Kennedy, B. J., Howard, C. J., and Chakoumakos, B. C. (1999).  Phase transitions in perovskite at elevated temperatures - a powder neutron diffraction study.  Journal of Physics: Condensed Matter, 11(6):1479–1488.
  * [Kieffer et al., 2020] Kieffer, J., Valls, V., Blanc, N., and Hennig, C. (2020).  New tools for calibrating diffraction setups.  Journal of Synchrotron Radiation, 27(2):558–566.
  * [Koch et al., 2021] Koch, R. J., Sinclair, R., McDonnell, M. T., Yu, R., Abeykoon, M., Tucker, M. G., Tsvelik, A. M., Billinge, S. J. L., Zhou, H. D., Yin, W.-G., and Bozin, E. S. (2021).  Dual Orbital Degeneracy Lifting in a Strongly Correlated Electron System.  Physical Review Letters, 126(18):186402.
  * [Kwei et al., 1995] Kwei, G. H., Billinge, S. J. L., Cheong, S.-W., and Saxton, J. G. (1995).  Pair-distribution functions of ferroelectric perovskites: Direct observation of structural ground states.  Ferroelectrics, 164(1):57–73.
  * [Kwei et al., 1993] Kwei, G. H., Lawson, A. C., Billinge, S. J. L., and Cheong, S. W. (1993).  Structures of the ferroelectric phases of barium titanate.  Journal of Physical Chemistry, 97(10):2368–2377.
  * [Müller et al., 1968] Müller, K. A., Berlinger, W., and Waldner, F. (1968).  Characteristic Structural Phase Transition in Perovskite-Type Compounds.  Physical Review Letters, 21(12):814–817.
  * [O’Keeffe and Hyde, 1977] O’Keeffe, M. and Hyde, B. G. (1977).  Some structures topologically related to cubic perovskite ($E2_{1}$), ReO3 ($D0_{9}$) and Cu3Au ($L1_{2}$).  Acta Crystallographica Section B, 33(12):3802–3813.
  * [Perez-Mato et al., 2010] Perez-Mato, J. M., Orobengoa, D., and Aroyo, M. I. (2010).  Mode crystallography of distorted structures.  Acta Crystallographica Section A, 66(5):558–590.
  * [Salje, 1990] Salje, E. (1990).  Phase transitions in ferroelastic and co-elastic crystals.  Ferroelectrics, 104(1):111–120.
  * [Sartbaeva et al., 2007] Sartbaeva, A., Wells, S. A., Thorpe, M. F., Bozin, E. S., and Billinge, S. J. L. (2007).  Quadrupolar Ordering in LaMnO3 Revealed from Scattering Data and Geometric Modeling.  Physical Review Letters, 99(15):155503.
  * [Sartbaeva et al., 2006] Sartbaeva, A., Wells, S. A., Thorpe, M. F., E. S. Bozin, and Billinge, S. J. L. (2006).  Geometric simulation of perovskite frameworks with Jahn-Teller distortions: application to the cubic manganites.  Physical Review Letters, 97:065501.
  * [Senn et al., 2016] Senn, M. S., Keen, D. A., Lucas, T. C. A., Hriljac, J. A., and Goodwin, A. L. (2016).  Emergence of Long-Range Order in BaTiO3 from Local Symmetry-Breaking Distortions.  Physical Review Letters, 116(20):207602.
  * [Skjærvø et al., 2019] Skjærvø, S. H., Meier, Q. N., Feygenson, M., Spaldin, N. A., Billinge, S. J. L., Bozin, E. S., and Selbach, S. M. (2019).  Unconventional Continuous Structural Disorder at the Order-Disorder Phase Transition in the Hexagonal Manganites.  Physical Review X, 9(3):031001.
  * [Wang et al., 2020] Wang, Z., Zhao, X.-G., Koch, R., Billinge, S. J. L., and Zunger, A. (2020).  Understanding electronic peculiarities in tetragonal FeSe as local structural symmetry breaking.  Physical Review B, 102(23):235121.
  * [Woodward, 1997a] Woodward, P. M. (1997a).  Octahedral Tilting in Perovskites. I. Geometrical Considerations.  Acta Cryst B, 53(1):32–43.
  * [Woodward, 1997b] Woodward, P. M. (1997b).  Octahedral Tilting in Perovskites. II. Structure Stabilizing Forces.  Acta Crystallographica Section B, 53(1):44–66.
  * [Yang et al., 2020] Yang, L., Koch, R. J., Zheng, H., Mitchell, J. F., Yin, W., Tucker, M. G., Billinge, S. J. L., and Bozin, E. S. (2020).  Two-orbital degeneracy lifted local precursor to a metal-insulator transition in MgTi2O4.  Physical Review B, 102(23):235128.
  * [Yashima and Ali, 2009] Yashima, M. and Ali, R. (2009).  Structural phase transition and octahedral tilting in the calcium titanate perovskite CaTiO3.  Solid State Ionics, 180(2):120–126.