# Representing Polymers as Periodic Graphs with Learned Descriptors for
Accurate Polymer Property Predictions

Evan R. Antoniuk  Peggy Li  Global Security Computing Applications Division,
Computing Directorate, Lawrence Livermore National Laboratory, Livermore, CA,
USA  Bhavya Kailkhura  Machine Intelligence Group/Center for Applied
Scientific Computing, Computing Directorate, Lawrence Livermore National
Laboratory, Livermore, CA, USA  Anna M. Hiszpanski

###### Abstract

One of the grand challenges of utilizing machine learning for the discovery of
innovative new polymers lies in the difficulty of accurately representing the
complex, multi-scale structures of polymeric materials. Although a wide array
of hand-designed polymer representations have been explored, there has yet to
be an ideal solution for the problems of how to capture the periodicity of
polymer structures, as well as how to automatically develop polymer
descriptors without the need for human feature design. In this work, we tackle
these problems through the development of our periodic polymer graph
representation. Specifically, our pipeline for polymer property predictions is
comprised of our polymer graph representation that naturally accounts for the
periodicity of polymers, followed by a message-passing neural network (MPNN)
that leverages the power of graph deep learning to automatically learn
chemically-relevant polymer descriptors. Across a diverse dataset of 10
polymer properties, we find that this polymer graph representation
consistently outperforms previous hand-designed representations with a sizable
20% average reduction in prediction error. Our results illustrate how the
incorporation of chemical intuition through directly encoding periodicity into
our polymer graph representation is responsible for a considerable improvement
in the accuracy and reliability of polymer property predictions. We also
demonstrate how combining polymer graph representations with message-passing
neural network architectures can automatically extract meaningful polymer
features that are consistent with human intuition, while achieving better
predictive performance than human-derived features. Whereas previous polymer
informatics have typically utilized hand-designed descriptors that have been
optimized for representing molecules, this work highlights the advancement in
predictive capability that is possible if using chemical descriptors that are
specifically optimized for capturing the unique chemical structure of
polymers.

## Introduction

The complex chemical structure of polymers provides chemists with an
unprecedented level of chemical design opportunities. As a result, efforts
focused on the rational design of polymers has led to their use in a
remarkably wide range of applications.[1, 2] By varying a polymer’s monomer
structure, its propagation site, or the polymer chain length, the properties
of the polymer can be tuned to suit the application of interest.[3, 4, 5, 6]
Formulations of polymers with chemical additives adds an even richer potential
design and functional space. However, efficiently navigating this vast
chemical design space to identify promising polymer materials for a given
application is incredibly challenging. Recently, the field of polymer
informatics has shown much promise in tackling this challenge by utilizing
machine learning (ML) models to learn the complex structure-property
relationships of polymers.[7, 8] Through the development of increasingly
accurate ML models, one can envision the future possibility of fully
computational polymer design workflows whereby polymers can be rapidly
designed in silico to satisfy multiple property criteria.

Accurate machine learning models require that the input to the model – be it
an image, text, or chemical structure – is transformed into a machine-readable
format that captures and represents all relevant features of the input. The
development of suitable polymer representations to use as inputs for machine
learning algorithms has been recognized as one of the grand challenges in the
field of polymer informatics.[9, 10] Unlike molecules that consist of a fixed
number of atoms and typically only contain dozens of atoms, a given polymer
can consist of a periodic chain of between 100-106 repeating chemical sub-
units.[11] This presents a considerable challenge for creating machine-
readable polymer representations since they must encode the periodic nature of
the polymer structure and also be compact enough to represent up to millions
of atoms. Despite these challenges, previous work on small molecules and
crystalline solids has demonstrated that developing physically-inspired
representations of atomistic systems is a fruitful direction towards realizing
significant improvements in property prediction accuracy.[12, 13, 14, 15]

In this work, we develop a new polymer representation that is the first to
explicitly capture the periodic nature of polymeric materials. First, we show
that the widely used monomer-based representations cannot sufficiently
describe the periodic, extended structure of polymers and can result in
unreliable property predictions. To solve this issue, we have developed a
novel polymer representation that we refer to as a polymer graph and that
naturally incorporates polymers’ periodicity in its construction. Across a
diverse spectrum of 10 unique polymer properties, we find that our polymer
graph representation consistently outperforms previous monomer-based
representations with a 20% average reduction in prediction error. Finally, we
demonstrate how our polymer graph representations can be combined with a
Message Passing Neural Network (MPNN) framework to automatically extract
chemically-relevant polymer descriptors directly from a dataset of polymer
properties. We find that these MPNN-learned polymer features yield more
accurate property predictions than predetermined features that have been hand-
crafted by chemists over years.

## Results

### Introduction to Polymer Graphs

Recently, a wide range of machine learning (ML) algorithms have been developed
to rapidly and accurately predict polymer properties in an a priori manner.
[16, 17] These previous approaches typically represent polymers by their
monomer repeat unit which is then converted into either i) molecular
descriptors that includes the presence of molecular fragments and/or property
descriptors (such as the molecular weight of the molecule)[18, 19, 20] or ii)
a molecular graph, where atoms are represented as nodes and chemical bonds are
represented as edges (Figure 1a)). [17] It should be noted that these
representations were first developed for molecules before being transferred to
the field of polymer informatics. Although these monomer-based representations
perform well at representing finite molecules, we show in this work that
monomer-based representations cannot adequately represent the extended chain
structures that define polymers.

Figure 1:  a) 5 polymer representations considered in this work. Monomer
graphs and chemical descriptors of the monomer have been reported previously,
whereas dimer graphs, trimer graphs, and polymer graphs are novel
representations developed in this work. b) Structure of the polymer graph
representation. All atoms (nodes) and bonds (edges) in the graph are
represented by learned feature vectors. In this example, the additional
periodic edge is featurized in the same manner as all other C-O edge vectors.
c) Representative graphs for two distinct polymers. Carbon and oxygen atoms
are depicted by black and red circles, respectively. For clarity, hydrogen
atoms are omitted. Generating monomer graphs does not produce a unique
representation and results in incorrect functional groups. Dimer and trimer
graph representations are unique for Polymers A and B, but result in incorrect
bonding environments on the atoms at the end of the chain. The polymer graph
representation developed in this work creates distinct representations for
Polymers A and B, while ensuring the correct bonding environment for terminal
atoms of the polymer chain and capturing the periodicity of the polymer.

Using monomer-based representations as inputs for ML algorithms results in at
least three notable issues: i) the generation of non-unique representations
for polymers that are composed of the similar monomer repeat units (see
"monomer graphs" of Polymer A and B in Figure 1 c), ii) incorrect bonding
environment of the polymer, and iii) a lack of periodicity. Polymers A and B
in Figure 1c) are two distinct polymers with notably different chemical
structures and properties. Since the monomer of a polymer is defined as the
smallest repeat unit of the polymer chain, the monomer of the two polymers in
Figure 1 c) is therefore the same C-C-O subunit (or CH3CH2OH if including
hydrogen atoms). Regardless of the specific details of the featurization used,
generating a representation of Polymer A and B based only on their monomers
will result in two identical representations (see Supplementary Note 1). In
turn, an ML algorithm that then uses these representations to learn a
structure-property relationship will, at best, only be able to learn the
average value of the property across the polymers, which may lead to large
prediction errors. The second issue that will arise from the use of monomer-
based representations is the incorrect featurization of the polymer. Polymer B
consists of a repeating (-C-C-O-) backbone, which is characterized by an ether
functional group (R-O-R, where R is any alkane). However, the monomer of
Polymer B ( C-C-O) does not contain an ether functional group since this ether
group is artificially fragmented when forming the monomer representation. This
is particularly problematic since both fingerprinting approaches and property
descriptors utilize the presence of functional groups and connectivity in
their representations (Supplementary Note 1).[21, 22, 23] Thirdly, the failure
to include periodicity in these monomer-based representations results in
representations that are inconsistent due to trivial translations along the
polymer backbone. The polymer (-C-C-O-)n can also be written as (-C-O-C-)n by
translating all atoms one position along the polymer chain. However, the
monomer representations for these two equivalent polymer representations would
be C-C-O and C-O-C, respectively. As shown in this example, the failure to
explicitly account for periodicity in these monomer-based representations can
result in inconsistent featurizations, which can negatively impact the ability
of ML algorithms to learn the underlying structure-property relationships.
These inconsistent featurizations are illustrated in Supplementary Note 1.

Fundamentally, all three of the aforementioned problems arise due to a failure
to account for the periodic nature of polymers. In other words, a finite
molecule (i.e. a monomer) cannot be used to represent a continuous, extended
chain of atoms without fragmenting chemical bonds in an unphysical manner. We
note that featurizing Polymers A and B by generating dimers (or any larger
oligomer chain) instead of monomers results in unique representations. Using a
dimer representation, however, still results in incorrect featurizations due
to the incomplete bonding environment of the atoms at the terminal end of the
dimer. In this work, we introduce a new polymer representation that we refer
to as a polymer graph that solves all three of these problems. The core idea
of the polymer graph is that adding a bond between the two terminal ends of
the polymer chain creates a periodic polymer representation (Figure 1b)). The
periodic nature of this representation eliminates the need for artificially
cleaving chemical bonds, ensuring that the correct bonding environment is
maintained for all atoms in the polymer. We note that this approach is
strongly reminiscent of the crystal graph representations that has had wide
success in representing periodic crystalline materials that extend
indefinitely in all three dimensions.[12] Crystal graph representations have
solved the problem of representing an infinite 3-D crystal by treating it as
the smallest repeat unit (referred to as the unit cell) that then transmits
chemical information by interacting with periodic copies of itself in all
directions. In an analogous fashion, the polymer graph developed in this work
is formed from the repeat unit of the polymer that is able to transmit
chemical information through the added periodic bond. We emphasize that this
periodic bond is used only for representing the polymer, but should not be
interpreted to mean that the polymers are always part of a chemical ring
system.

Graph-based representations of molecules and crystalline solids have recently
been the subject of considerable research efforts and are currently the state-
of-the-art representation for predicting the properties of these systems.[24,
25, 13] Generally these graphs represent each atom in the chemical system,
$i$, by an atom embedding vector (hi) which is updated through messages that
are sent through the edges of the graph. These update messages typically
operate on all atoms and bonds in the local neighborhood of atom $i$, where
the local neighborhood is given by the edges of the graph that are connected
to atom $i$. To generate our polymer graph representation, we first create a
standard trimer graph and then add an additional edge to connect the two
terminal atoms of the trimer (Figure 1 b)). This added edge in our polymer
graph representation not only maintains the correct bonding environment for
all atoms, but also serves to ensure that messages (i.e. chemical information)
can be sent throughout the polymer graph. For example, the left-most carbon
atom in the monomer graph in Figure 1a) will only directly receive messages
from the central carbon atom since this is the only atom that it is bonded to.
This is not an accurate representation of the true chemical structure of this
polymer since this carbon atom is directly bonded to both a carbon atom and an
oxygen atom. On the other hand, since our polymer graph representation creates
an edge between the terminal carbon and oxygen atoms, messages are directly
passed between both of these two atoms, thus allowing the atom embedding
vector of the terminal carbon atom to more accurately represent its true
chemical environment. We construct the polymer graph from a trimer graph
(rather than a dimer or monomer graph) to ensure that the added periodic bond
is between atoms that are distinct and not already bonded to each other.

### Performance of Graph Representations

To demonstrate the accuracy and generality of our polymer graph representation
for predicting polymer properties, we collect a dataset of polymer properties
containing 15,219 datapoints across 10 distinct properties (see Methods for
dataset details). Notably, this dataset includes a diverse set of property
types (thermal, structural, optical, electronic) as well as data collected
from both experimental and computational sources, which allows us to assess
the capability of our approach to learn meaningful polymer representations for
general property prediction tasks. To perform a baseline comparison, we also
train ML models on all 10 properties with two commonly used monomer-based
representations: i) molecular descriptors and ii) molecular graphs of the
monomer (Table 1). All graph-based models are trained with a MPNN model,
whereas a random forest (RF) model is used for training with the monomer
descriptors (see Methods for model architecture and hyperparameter details).

Table 1: Performance of random forest (RF) and message-passing neural network (MPNN) models for predicting polymer properties. For each property, the reported error is the root mean square error (RMSE) on the validation set of 5-fold cross validation. The mean normalized RMSE is obtained by normalizing the RMSE of a model against the RMSE of the RF model and then averaging across all 10 properties. The best performing model for each property is shown in bold. Property |  Unit |  Number of Data Points |  RF/ Molecular Descriptors |  MPNN/ Monomer Graph |  MPNN/ Polymer Graph
---|---|---|---|---|---
Atomization Energy | eV/atom | 390 | 0.100 | 0.057 | 0.038
Crystallization Tendency | % | 432 | 18.050 | 19.850 | 20.180
Band Gap Chain | eV | 3380 | 0.579 | 0.528 | 0.449
Band Gap Bulk | eV | 561 | 0.698 | 0.581 | 0.528
Electron Affinity | eV | 368 | 0.450 | 0.306 | 0.289
Ionization Energy | eV | 370 | 0.456 | 0.473 | 0.385
Refractive Index | - | 382 | 0.104 | 0.091 | 0.090
Dielectric Constant | - | 382 | 0.593 | 0.533 | 0.535
Glass Transition Temperature | K | 7235 | 41.39 | 38.79 | 34.91
Density | g/cc | 1719 | 0.099 | 0.090 | 0.086
Mean Normalized RMSE | - | - | 1.000 | 0.875 | 0.800

Comparing the performance of the RF/Molecular Descriptors and the MPNN/Monomer
Graph in Table 1 illustrates the increase in performance that is realized when
using graph-based models for polymer property prediction. Across 8 out of 10
properties, the MPNN/Monomer graph model outperforms the RF/Molecular
Descriptors model. On average, the MPNN/Monomer Graph model achieves 12 %
lower error than the RF/Molecular Descriptors model, with remarkably large
reductions in error for predicting atomization energy (43%), band gap bulk
(17%), and electron affinity (32%). Both models perform similarly at
predicting ionization energy, whereas the RF/Molecular Descriptors model
performs slightly better only at predicting the crystallization tendency. The
poor performance of the monomer graph model for predicting crystallization
tendency may be due to the fact that these values were calculated with the use
of an approximate group contribution method.[26] The low fidelity of this
dataset may therefore be preventing the graph model from learning reasonable
atom/bond features for predicting crystallization tendency. Based on the
results across all properties, we conclude that graph-based models can
generally be expected to perform, at worst, comparably to RF/Molecular
Descriptor models, with the potential to perform significantly better on
certain properties. This finding is in agreement with a previous report which
showed that monomer graph-based models outperform ECFP molecular descriptors
for predicting the glass transition temperature, melting temperature, density,
and elastic modulus for a subclass of polymers known as polyamides. [17] In
the present work, we expand this finding to show that graph-based models can
generally be expected to outperform models based on molecular descriptors
across a broader range of polymer properties and for datasets that contain all
subclasses of polymers.

### Role of Graph Construction on Predictions

Now that we have demonstrated the ability of graph-based models to perform
well at predicting polymer properties, we now explore how the construction of
the graph itself can influence model performance. Table 1 illustrates the
performance of monomer graph and polymer graph representations with the same
MPNN model. Although we have already demonstrated that the MPNN/Monomer Graph
model considerably outperforms the current state-of-the-art RF/Molecular
Descriptor models, we find that using our polymer graph representation results
in a further 10% average reduction in error compared to the monomer graph
model and does not perform significantly worse than the monomer graph model
for any of the properties. The polymer graph representation results in a
significant reduction in error for predicting atomization error (33%), band
gap chain (15%), ionization energy (19%) and glass transition temperature
(10%), compared to the monomer graph representation. Since the improvement in
property prediction accuracy achieved by the polymer graph is highly property-
dependant, we note that the 10% average reduction in error reported here is
likely to vary depending on the chosen set of properties of interest.
Nevertheless, we do not find any properties where the polymer graph performs
significantly worse than the monomer graph representation, leading us to
conclude that the polymer graph should be expected to yield more accurate
property predictions across a wide range of properties.

Our original hypothesis to motivate the development of our polymer graph
representation was that it can improve over monomer graph representations by
being able to distinguish between polymers that have the same monomer
representation. To test this idea, we construct a test set of 179 glass
transition temperature (Tg) datapoints where each polymer in the dataset has
at least one other polymer that contains the same monomer representation, as
depicted by polymers A-C in Figure 2. We train a MPNN/Monomer Graph model and
a MPNN/Polymer Graph model on a Tg training set that does not include these
179 polymers and then evaluate their ability to predict the Tg of the 179
polymers in the test set. This experiment therefore allows us to quantify the
impact on prediction accuracy that results from the inability of monomer
graphs to distinguish between polymers that have the same monomer
representation.

Figure 2: Comparison of the performance of monomer graph and polymer graph
representations for predicting Tg of a test set of polymers that contain at
least one other polymer with the same monomer structure.

As illustrated in Figure 2, the inability of the MPNN/Monomer graph model to
distinguish between polymers A-C results in the same Tg prediction for all
three polymers. Since the experimental Tg value varies largely between these
three polymers, this results in an unacceptably large prediction error of over
60∘C. Although polymers A-C are chosen as an illustrative example, a more
statistically significant comparison can be obtained from the performance on
the entire test set of 179 polymers. Consistent with our original hypothesis,
we see that the inability of the monomer graph representation to distinguish
between the polymers in this test set results in significantly worse
performance when compared to Table 1 (49∘C vs. 39∘C). On the other hand, the
polymer graph representation is not only able to distinguish between these
polymers, but also quantitatively capture the trends in Tg values for polymers
A-C with a RMSE that is nearly 60% less than that of the monomer graph. These
results therefore provide a partial explanation for the superior of
performance of the polymer graph representation in Table 1, relative to the
monomer graph representation. By being able to distinguish between polymers of
the same monomer structure, the polymer graph representation is therefore
better able to learn quantitative structure-property relationships for the
properties in Table 1. Since the proportion of polymers that have the same
monomer structure is relatively small ( 2% of the Tg dataset), we anticipate
that this effect is currently small, but will grow in importance as more
polymers are discovered.

The results of the experiment shown in Figure 2 also have important practical
implications for the reliability that can be expected when using ML methods to
discover new polymers. The Tg of a polymer is the temperature at which the
thermal motions of the polymers are large enough to overcome the forces that
hold the amorphous polymer together, such that the polymer is able to undergo
large-scale molecular motions.[27] Accurately predicting the glass transition
temperature is therefore of great importance since Tg is often treated as the
maximum temperature that the polymer can be used.[28] Despite having identical
chemical compositions, the three polymers shown in Figure 2 exhibit a
remarkably broad range of experimental Tg values due to their unique chemical
structures. In particular, the presence of the hydroxyl (OH) group on Polymer
A leads to an increase in interchain attractive forces, resulting in a
significantly larger Tg value than for Polymers B and C.[29] The inability of
monomer-based graph representations to distinguish between these polymers
therefore severely hurts the trustworthiness of the model. If one were to use
the MPNN/Monomer Graph model to discover new polymers, the large errors in
this subset of Tg predictions could manifest in the discovery of polymers that
are thermodynamically unstable at the desired operating temperature. By
comparison, the consistent performance seen by the MPNN/Polymer Graph model
helps to instill confidence that the predictions will be robust across a large
dataset of polymer materials.

### Ablation Studies on Periodic Graphs

Although the results in Table 1 and Figure 2 suggest that the polymer graph
representation can more accurately predict polymer properties than both the
monomer graph and molecular descriptor representations, they do not fully
explain what is causing this improvement in performance. Understanding
specifically how polymer representation can impact predictive performance is
important for the development of increasingly better polymer representations.
To accomplish this goal, we perform several ablation experiments whereby we
systematically alter components of the polymer graph representation to allow
us to quantify the contribution of each component to overall model performance
(Figure 3). We conduct this experiment by training the same MPNN model as
before on the same 10 properties listed in Table 1. For clarity, we illustrate
the performance of each representation by calculating the average performance
across all properties (Figure 3). The property-specific performance of each
representation is given in Table S1.

Figure 3: Results of ablation experiments for various graph-based polymer
representations. The mean normalized RMSE is calculated for all 10 properties
listed in Table 1 and normalized against the performance of the RF/Molecular
Descriptors representation. All property-specific RMSE values are given in
Table S1.

The first ablation experiment that we perform is to remove the periodic edge
in the polymer graph representation, which results in a trimer graph. This
ablation experiment highlights the importance of including periodicity in
polymer representations. Although the polymer graph and trimer graph
representations only differ by a single edge in their graphs, the inclusion of
periodicity in the polymer graph representation results in a 4% average
reduction in performance error. Most notably, a 32% reduction in the error of
atomization energy predictions is observed when using the polymer graph
representation instead of the trimer graph (Table S1). Next, we explore the
effect of increasing the size of the polymer chain that is used in the graphs
by comparing the performance achieved when using the trimer graph and dimer
graph representations. We do not observe any significant performance
difference between these two representations. This result is perhaps
unsurprising since the trimer graph does not provide any additional chemical
information compared to the dimer graph representation. Based on this finding,
we anticipate that using even larger graph representations (i.e. tetramer or
pentamer) will not improve model performance, but will unnecessarily require
more computational resources for model training. Finally, comparing the
performance of the dimer graph and the monomer graph representations
illustrates that including more than one repeat unit in the graph
significantly improves performance. It may seem counter-intuitive that the
trimer graph and dimer graph representations perform similarly, whereas dimer
graph representations considerably outperform monomer graph representations.
However, this result can be explained by noting that the dimer graph
representation contains additional chemical bonding information, relative to
the monomer graph, since it directly represents the atoms that are bonded
together to form the polymer chain. On the other hand, the trimer/dimer graph
representations do not differ in the chemical information that they contain,
but are simply scaled versions of each other.

### Analysis of MPNN Learned Polymer Descriptors, Hand-Designed Descriptors
and Chemical Intuition

We note from the ablation experiments illustrated in Figure 3 that the largest
reduction in error relative to the RF/Monomer Descriptors baseline is due to
using graph-based representations with the MPNN model. Using graph-based
models results in at least a 13% average reduction in error, relative to the
RF/Monomer Descriptors baseline, whereas further refinements to the graph
representation only affects the average performance by approximately 3-5%. To
better understand why using graph-based representations in combination with a
MPNN architecture leads to such a significant performance improvement, we
analyze the learned feature vectors that are generated from the MPNN model.
The MPNN model utilized in this work uses learned bond and atom features that
are aggregated into a single fixed-length vector. This aggregated vector is
then passed through a feed-forward neural network with two hidden layers to
generate the final property prediction. As a result, the vector in the final
hidden layer of this feed-forward neural network can be considered as a
learned encoding of the polymer that has been optimized by the MPNN model to
make accurate property predictions. These vector encodings can therefore be
analyzed to determine what chemical information the MPNN model is using to
generate its predictions.

Whereas the RF/Monomer Descriptors baseline used previously (Table 1) utilizes
pre-determined molecular descriptors (such as molecular weight), the vector
encodings extracted from our MPNN model are vector representations that have
been specifically optimized for the particular property that they are trained
on. Comparing the performance of the RF/Monomer Descriptors baseline and a RF
model trained with these learned encodings therefore allows us to
quantitatively understand if the superior performance of the graph-based
models are due to these optimized encodings or the MPNN model architecture. To
perform this experiment, we train the MPNN model on k-1 folds (k=5) of the Tg
dataset. This trained MPNN model is then used to generate vector encodings for
the polymers in the same k-1 folds of the dataset. Then, these vector
encodings are used as inputs to train the RF model on the same k-1 folds of
the dataset that the MPNN model was trained on. Finally, the test performance
of the RF model is calculated for the remaining fold of the dataset that was
not seen by the MPNN model nor the RF model in training. This procedure is
repeated for k iterations to ensure that each fold of the dataset is used as
the test set once and the reported performance is averaged across these k
iterations. This training procedure prevents data leakage by ensuring that the
test set is not seen by the MPNN model or the RF model before determining test
performance.

The results of this experiment are given in Figure 4a). Remarkably, we find
that the RF model with learned descriptors performs comparably to the
MPNN/Polymer Graph model at predicting Tg. This result suggests that the
excellent performance of the MPNN/Polymer Graph model is not due to MPNNs
being a more powerful model architecture. Rather, the learned descriptors that
the MPNN model extracts provides a more meaningful representation of polymers,
leading to superior predictive capability. This finding is exciting because it
suggests that the use of learned polymer representations that better capture
the underlying chemical properties of polymers can lead to greatly improved
predictive capabilities, even with relatively simple model architectures like
random forests. It is important to note that the predetermined molecular
descriptors that have been previously used for polymer prediction tasks were
first developed as general descriptors for molecules, but with no guarantee
that these descriptors will contain properties that are relevant for
describing polymer properties.[21] On the other hand, these learned MPNN
descriptors are specifically trained to capture the relevant components of a
polymer’s structure that give rise to its properties. In this way, the MPNN
model explored in this work provides a recipe for automatically extracting
meaningful polymer representations directly from polymer property datasets in
a manner that does not require any hand-crafted featurization.

Although we have shown that the MPNN model is capable of extracting features
that achieve better performance than predetermined features, it is nonetheless
important to gain an understanding of what the MPNN learned features
physically represent. To gain this chemical insight, we calculate the
correlation between the 5 heaviest weighted features of the final layer of the
MPNN/Polymer Graph model and the predetermined molecular descriptors for the
polymers in the Tg dataset. The 5 molecular descriptors that are most
correlated with the MPNN learned features for the MPNN/Polymer Graph model are
depicted in Figure 4b). Two of the most important known structural features
for predicting Tg in polymers is: i) the presence of polar functional groups
such as NH2 and OH that increase intermolecular forces, resulting in a higher
Tg value and ii) the stiffness of the polymer chain, which is typically
described in terms of the presence of rigid ring structures in the polymer
backbone.[27] The molecular descriptors shown in Figure 4b) are colored
according to which one of these two known structural factors they explicitly
capture. Interestingly, we see that the majority of the learned MPNN features
are strongly correlated with these known relevant chemical features, which
suggests that our MPNN model is autonomously learning the chemical features
that are responsible for Tg directly from the dataset of Tg values. In
particular, the most heavily weighted MPNN feature is correlated with the
presence of polar groups, whereas the other top 5 features provide complex
descriptions of the presence of ring structures in the polymer. However, we
emphasize that these MPNN polymer encodings also contain new information that
is not captured by the predetermined descriptors, as is shown by the superior
performance of the RF model when using learned descriptors instead of
predetermined descriptors (Figure 4a).

Figure 4:  a) Comparison of RF models with predetermined molecular descriptors
and learned descriptors against the MPNN/Polymer Graph model for predicting
glass transition temperature. b) Correlation between the 5 most heavily
weighted features in the last layer of the MPNN model and the predetermined
molecular descriptors. Features are colored according to the type of chemical
information that they capture. For clarity, the molecular descriptors from
Rdkit have been renamed to more explicitly describe their property. The
original Rdkit molecular descriptors names and an explanation of these
descriptors are given in Supplementary Note 2.

## Discussion

Generally, this work has demonstrated how the development of representations
that better account for the unique and complex chemical structure of polymers
can be a powerful strategy for increasing the accuracy of polymer property
predictions. Through a number of illustrative examples, we have demonstrated
how a failure to account for periodicity in polymers can lead to property
predictions that are inaccurate and unreliable. We have shown that these
issues can be mitigated by our novel polymer graph representation that
explicitly accounts for the inherent periodicity of polymeric materials. The
results presented here therefore demonstrate that utilizing chemical intuition
to construct physically-informed polymer representations is a viable strategy
for advancing the predictive capabilities of ML models.

The polymer informatics community has benefited greatly from the adaptation of
machine learning algorithms and representations that were first developed for
molecules. This approach is based on the assumption that general molecular
descriptors are likely to work well across a broad spectrum of chemical
classes, including polymers. In this work, we test this assumption by
exploring the improvement in property predictions that can be achieved through
the use of learned descriptors that are specifically optimized to capture the
unique structural characteristics of polymers. We have shown that utilizing
polymer-specific descriptors results in considerable improvements in the
accuracy of polymer property predictions, compared to standard molecular
descriptors. These results suggest that the structural characteristics that
give rise to polymer properties are distinct from general molecular
descriptors. This is perhaps unsurprising since it is well known that polymers
are affected by structural elements that do not have molecular analogues, such
as the branching structure of the polymer or its interchain interactions. We
anticipate that the continued development of representations that are
customized to capture the unique and complex structure of polymeric systems is
therefore a fruitful direction for further pushing the limits of predicting
polymer properties with ML.

## Methods

### Data

The polymer dataset used in this work consists of 2 experimentally measured
properties that were manually obtained from the PolyInfo database (glass
transition temperature (Tg), and density.[30] The data for the other 8
properties (atomization energy, crystallization tendency, band gap chain, band
gap bulk, electron affinity, ionization energy, refractive index, dielectric
constant) were previously calculated with density functional theory (DFT).[18]
Altogether, this dataset consists of 15,219 total datapoints, spanning 9,935
unique homopolymer structures. For the experimentally measured properties, we
exclude all measurements that correspond to compound or composite polymer
samples. After these measurements are excluded, the experimental property
value for each polymer is taken to be the median of all measurements.

### Molecular Descriptors

The molecular descriptors used in combination with the random forest model
were obtained from the open source cheminformatics tookit, RDKit.[21] In
total, these descriptors contain 200 unique features: 115 of which describe
chemical characteristics (such as the number of rings or molecular weight),
while the other 85 are binary features that describe the presence of molecular
fragments (such as the presence of an amide functional group). For each
property, any features that have a constant value across the whole dataset are
removed prior to training. The molecular descriptors are generated from the
monomer SMILES string of the polymer that has the ’*’ characters removed. For
example, the monomer input SMILES string for the polymer (-CH2-CH2-O-)n would
be ’CCO’.

### Models

The Message Passing Neural Network (MPNN) used in this work is implemented in
the Chemprop package.[13] Notably, this network utilizes directed bond-level
messages (rather than the atom-level messages used in traditional MPNNs) to
transmit chemical information across the molecule. This network architecture
operates on input molecular graphs to first generate a vector encoding of the
graph, followed by a feed-forward neural network that converts the graph
encodings into property predictions. The size of the bond message vectors is
set to 300 and a total of 3 message-passing steps are used. The feed-forward
neural network consists of 2 hidden layers, each consisting of 300 hidden
neurons. All hyperparameters used are the Chemprop default values. We point
interested readers to the Chemprop documentation for an extensive list of all
hyperparameters used in this architecture.[13] Unless otherwise stated, all
models are trained with 5-fold cross-validation (CV) for 50 epochs on each
fold. The best-performing model is chosen by selecting the model with the
lowest root-mean-square error (RMSE) performance on the validation set. The
Random Forest (RF) model was implemented in scikit-learn version 0.24.1 with
500 trees.[31]

## References

  * [1] Hou, W., Xiao, Y., Han, G. & Lin, J.-Y.  The Applications of Polymers in Solar Cells: A Review.  _Polymers_ 11, 143, DOI: 10.3390/polym11010143 (2019).  Number: 1 Publisher: Multidisciplinary Digital Publishing Institute.
  * [2] Briscoe, B. J. & Sinha, S. K.  CHAPTER 1 - Tribological applications of polymers and their composites: Past, present and future prospects.  In Friedrich, K. & Schlarb, A. K. (eds.) _Tribology and Interface Engineering Series_ , vol. 55 of _Tribology of Polymeric Nanocomposites_ , 1–14, DOI: 10.1016/S1572-3364(08)55001-4 (Elsevier, 2008).
  * [3] Sharma, V. _et al._ Rational design of all organic polymer dielectrics.  _Nature Communications_ 5, 4845, DOI: 10.1038/ncomms5845 (2014).  Number: 1 Publisher: Nature Publishing Group.
  * [4] Pacułt, J., Rams-Baron, M., Chrząszcz, B., Jachowicz, R. & Paluch, M.  Effect of Polymer Chain Length on the Physical Stability of Amorphous Drug–Polymer Blends at Ambient Pressure.  _Molecular Pharmaceutics_ 15, 2807–2815, DOI: 10.1021/acs.molpharmaceut.8b00312 (2018).  Publisher: American Chemical Society.
  * [5] Mun, J. _et al._ Conjugated Carbon Cyclic Nanorings as Additives for Intrinsically Stretchable Semiconducting Polymers.  _Advanced Materials_ 31, 1903912, DOI: 10.1002/adma.201903912 (2019).  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/adma.201903912.
  * [6] Selivanova, M. _et al._ Branched Polyethylene as a Plasticizing Additive to Modulate the Mechanical Properties of pi-Conjugated Polymers.  _Macromolecules_ 52, 7870–7877, DOI: 10.1021/acs.macromol.9b01697 (2019).  Publisher: American Chemical Society.
  * [7] Goswami, S., Ghosh, R., Neog, A. & Das, B.  Deep learning based approach for prediction of glass transition temperature in polymers.  _Materials Today: Proceedings_ 46, 5838–5843, DOI: 10.1016/j.matpr.2021.02.730 (2021).
  * [8] Tao, L., Varshney, V. & Li, Y.  Benchmarking Machine Learning Models for Polymer Informatics: An Example of Glass Transition Temperature.  _Journal of Chemical Information and Modeling_ 61, 5395–5413, DOI: 10.1021/acs.jcim.1c01031 (2021).  Publisher: American Chemical Society.
  * [9] Chen, L. _et al._ Polymer informatics: Current status and critical next steps.  _Materials Science and Engineering: R: Reports_ 144, 100595, DOI: 10.1016/j.mser.2020.100595 (2021).
  * [10] Patra, T. K.  Data-Driven Methods for Accelerating Polymer Design.  _ACS Polymers Au_ 2, 8–26, DOI: 10.1021/acspolymersau.1c00035 (2022).  Publisher: American Chemical Society.
  * [11] Yamauchi, Y. _et al._ Two-Step Divergent Synthesis of Monodisperse and Ultra-Long Bottlebrush Polymers from an Easily Purifiable ROMP Monomer.  _Angewandte Chemie International Edition_ 60, 1528–1534, DOI: 10.1002/anie.202009759 (2021).  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/anie.202009759.
  * [12] Xie, T. & Grossman, J. C.  Crystal Graph Convolutional Neural Networks for an Accurate and Interpretable Prediction of Material Properties.  _Physical Review Letters_ 120, 145301, DOI: 10.1103/PhysRevLett.120.145301 (2018).  Publisher: American Physical Society.
  * [13] Yang, K. _et al._ Analyzing Learned Molecular Representations for Property Prediction.  _Journal of Chemical Information and Modeling_ 59, 3370–3388, DOI: 10.1021/acs.jcim.9b00237 (2019).  Publisher: American Chemical Society.
  * [14] Schindler, P., Antoniuk, E. R., Cheon, G., Zhu, Y. & Reed, E. J.  Discovery of materials with extreme work functions by high-throughput density functional theory and machine learning.  _arXiv:2011.10905 [cond-mat, physics:physics]_ (2020).  ArXiv: 2011.10905.
  * [15] Nguyen, P. _et al._ Predicting Energetics Materials’ Crystalline Density from Chemical Structure by Machine Learning.  _Journal of Chemical Information and Modeling_ 61, 2147–2158, DOI: 10.1021/acs.jcim.0c01318 (2021).  Publisher: American Chemical Society.
  * [16] Huan, T. D. _et al._ A polymer dataset for accelerated property prediction and design.  _Scientific Data_ 3, 160012, DOI: 10.1038/sdata.2016.12 (2016).  Number: 1 Publisher: Nature Publishing Group.
  * [17] Park, J. _et al._ Prediction and Interpretation of Polymer Properties Using the Graph Convolutional Network.  _ACS Polymers Au_ DOI: 10.1021/acspolymersau.1c00050 (2022).  Publisher: American Chemical Society.
  * [18] Kuenneth, C. _et al._ Polymer informatics with multi-task learning.  _Patterns_ 2, 100238, DOI: https://doi.org/10.1016/j.patter.2021.100238 (2021).
  * [19] Wu, S. _et al._ Machine-learning-assisted discovery of polymers with high thermal conductivity using a molecular design algorithm.  _npj Computational Materials_ 5, 1–11, DOI: 10.1038/s41524-019-0203-2 (2019).  Number: 1 Publisher: Nature Publishing Group.
  * [20] Barnett, J. W. _et al._ Designing exceptional gas-separation polymer membranes using machine learning.  _Science Advances_ 6, eaaz4301, DOI: 10.1126/sciadv.aaz4301 (2020).  Publisher: American Association for the Advancement of Science.
  * [21] Landrum, G.  RDKit: Open-Source Cheminformatics.  http://www.rdkit.org (2006).
  * [22] Morgan, H. L.  The Generation of a Unique Machine Description for Chemical Structures-A Technique Developed at Chemical Abstracts Service.  _Journal of Chemical Documentation_ 5, 107–113, DOI: 10.1021/c160017a018 (1965).  Publisher: American Chemical Society.
  * [23] Rogers, D. & Hahn, M.  Extended-Connectivity Fingerprints.  _Journal of Chemical Information and Modeling_ 50, 742–754, DOI: 10.1021/ci100050t (2010).  Publisher: American Chemical Society.
  * [24] Gasteiger, J., Becker, F. & Günnemann, S.  GemNet: Universal Directional Graph Neural Networks for Molecules.  _arXiv:2106.08903 [physics, stat]_ (2022).  ArXiv: 2106.08903.
  * [25] Gasteiger, J., Groß, J. & Günnemann, S.  Directional Message Passing for Molecular Graphs.  _arXiv:2003.03123 [physics, stat]_ (2022).  ArXiv: 2003.03123.
  * [26] Venkatram, S. _et al._ Predicting Crystallization Tendency of Polymers Using Multifidelity Information Fusion and Machine Learning.  _The Journal of Physical Chemistry B_ 124, 6046–6054, DOI: 10.1021/acs.jpcb.0c01865 (2020).  Publisher: American Chemical Society.
  * [27] Bicerano, J.  _Prediction of Polymer Properties_ (CRC Press, 2002).  Google-Books-ID: 5L1USsiq6KcC.
  * [28] Hale, A., Macosko, C. W. & Bair, H. E.  Glass transition temperature as a function of conversion in thermosetting polymers.  _Macromolecules_ 24, 2610–2621, DOI: 10.1021/ma00009a072 (1991).  Publisher: American Chemical Society.
  * [29] Feldstein, M. M., Shandryuk, G. A. & Platé, N. A.  Relation of glass transition temperature to the hydrogen-bonding degree and energy in poly(N-vinyl pyrrolidone) blends with hydroxyl-containing plasticizers. Part 1. Effects of hydroxyl group number in plasticizer molecule.  _Polymer_ 42, 971–979, DOI: 10.1016/S0032-3861(00)00445-6 (2001).
  * [30] Otsuka, S., Kuwajima, I., Hosoya, J., Xu, Y. & Yamazaki, M.  PoLyInfo: Polymer Database for Polymeric Materials Design.  _2011 International Conference on Emerging Intelligent Data and Web Technologies_ DOI: 10.1109/EIDWT.2011.13 (2011).
  * [31] Pedregosa, F. _et al._ Scikit-learn: Machine learning in Python.  _Journal of Machine Learning Research_ 12, 2825–2830 (2011).

## Acknowledgements

This work was performed under the auspices of the U.S. Department of Energy by
Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.

## Supplementary Information

### Supplementary Note 1: Featurization of Monomers

Figure S1. Visualization of molecular descriptors generated in RDKit for 4
representative SMILES strings for polyethylene glycol, (-CH2-CH2-O-)n. The
molecular property features are normalized by the maximum value that they have
in the full dataset of Tg values.

To illustrate the inconsistencies that arise when using preexisting molecular
descriptors to represent polymers, Figure S1 provides a visualize of the
molecular descriptors generated from RDKit for 4 representative SMILES
strings. According to the syntax of the SMILES strings, the ’*’ character is
used to indicate the connection between repeating units of a polymer. The
SMILES notation also does not explicitly indicate hydrogen atoms, but instead
assumes that all atoms have the number of hydrogen atoms required for full
valency. The ’*CCO*’ SMILES string therefore represents a (-CH2-CH2-O-)n
polymer and the ’*COC*’ SMILES string represents the same polymer,
(-CH2-O-CH2-)n, that has been trivially translated along the polymer backbone.
Even though ’*CCO*’ and ’*COC*’ represent the same polymer, Figure S1
illustrates the different feature vectors that are generated from these SMILES
strings. An understanding of why this difference arises can be obtained by
looking at the molecular fragment features that are generated for these two
SMILES string. The only difference in the molecular fragment features is the
’fr_ether’ feature, which is a binary feature that indicates the presence of
an ether functional group. As shown in Figure S1, the ’*CCO*’ has fr_ether=0
(indicating the absence of an ether), whereas ’*COC*’ has fr_ether=1
(indicating the presence of an ether). Evidently, the inconsistent
featurization between these two chemically identical SMILES strings is due to
a lack of periodicity. Since the ’*CCO*’ SMILES string does not treat the
terminal oxygen atom as being bonded to two neighboring carbon atoms, using
the ’*CCO*’ SMILES string results in an incorrect featurization of fr_ether=0,
when this polymer does in fact contain an ether.

The impact of including the ’*’ character in the SMILES representation is
illustrated with feature vectors of the ’CCO’ and ’COC’ SMILES strings, which
are the monomer subunits of ’*CCO*’ and ’*COC*’, respectively. Interestingly,
all four of these inputs yield four different feature vectors. Although these
inconsistencies could be accounted for by designating a standardized SMILES
string to use as an input for a given polymer, it is unlikely that such a
standardized system would always capture the correct periodicity of the
polymer. Instead, the use of the polymer graph representation as described in
the main text ensures that the periodicity of the polymer is maintained and
always generates a consistent representation.

### Supplementary Note 2: Explanation of Features in Figure 4

Figure S2. Same as Figure 4b, but molecular features are named according to
the Rdkit documentation.

Here, we list the original names of the features (as described in the RDKit
documentation). Each feature in the main text is also categorized into
capturing 1) the presence of polar groups, 2) the number of rings, and 3)
other. For the features that do not explicitly represent these groups, we also
provide a simple justification for classifying each feature into each group.

fr_NH1: The presence of a NH1 group in the molecule. Polar Groups.

NumHDonors: Number of hydrogen bond donors in the molecule. Polar Groups.

NHOHCount: The total number of NH or OH functional groups in the molecule.
Polar Groups.

Fr_aniline: The presence of an aniline functional group in the molecule. Since
aniline contains both a polar group (NH) as well as a ring, we classify this
feature as capturing both 1) presence of polar groups and 2) the number of
rings.

VSA_EState2: The electrotopological state of an atom, summed for all atoms in
the molecule for a van der Waals surface area between 4.78-5.00. Other.

RingCount: Number of rings in the molecule. Number of rings.

NumAromaticRings: Number of aromatic rings in the molecule. Number of rings.

BertzCT: A topological descriptor that describes the ’complexity’ of a
molecule. Other.

PEOE_VSA14: Approximate descriptor for describing a molecules surface area.
Other.

Fr_benzene: The presence of a benzene ring in the molecule. Number of rings.

NumAromaticCarbocycles: Number of aromatic carbocyles in the molecule. Number
of rings.

VSA_EState6: Same as VSA_EState2, but calculated for a van der Waals surface
area between 6.00-6.07. Other.

SlogP_VSA8: The log P value of a molecule is a partition coefficient that
describes its lipophilicity. This descriptor describes the number of atoms in
the surface area of molecule that have log P values between 0.25 and 0.30. We
include this feature in Group 1 since the lipophilicity of a molecule is
heavily impacted by the polarity of a molecule.

Fr_bicyclic: Presence of bicyclic substructures in the molecule. Number of
rings.

FractionCSP3: Fraction of carbons in the molecule that are SP3 hybridized.
Other.

SMR_VSA7: Molar refractivity (MR) describes the polarizability of a molecule.
Similar to the SlogP_VSA8 feature, this descriptor describes the number of
atoms in the surface area of a molecule that have estimated molar refractivity
values between 3.05-3.63. Other.

Table S1: Performance of polymer graph, trimer graph, dimer graph and monomer graph representations for all 10 polymer properties. All representations are trained with the MPNN architecture described in the Methods section. Property | Unit | Num. Data Points | Polymer Graph | Trimer Graph | Dimer Graph | Monomer Graph
---|---|---|---|---|---|---
Atomization Energy | eV/atom | 390 | 0.038 | 0.056 | 0.055 | 0.057
Crystallization Tendency | % | 432 | 20.18 | 21.14 | 21.06 | 19.85
Band Gap Chain | eV | 3380 | 0.449 | 0.450 | 0.453 | 0.528
Band Gap Bulk | eV | 561 | 0.528 | 0.527 | 0.532 | 0.581
Electron Affinity | eV | 368 | 0.289 | 0.275 | 0.282 | 0.306
Ionization Energy | eV | 370 | 0.385 | 0.435 | 0.435 | 0.473
Refractive Index | - | 382 | 0.090 | 0.088 | 0.088 | 0.091
Dielectric Constant | - | 382 | 0.535 | 0.520 | 0.529 | 0.533
Glass Transition Temperature | K | 7235 | 34.91 | 37.00 | 34.79 | 38.79
Density | g/cc | 1719 | 0.086 | 0.088 | 0.086 | 0.090


# Microscopic Origin Of Room Temperature Ferromagnetism in a Double Perovskite
Sr2FeReO6: a first principle and model Hamiltonian study

Shishir K. Pandey  Department of Condensed Matter Physics and Material
Science, S. N. Bose National Center for Basic Sciences,Kolkata-700106  Ashis
K. Nandy  Department of Condensed Matter Physics and Material Science, S. N.
Bose National Center for Basic Sciences,Kolkata-700106  School of Physical
Sciences, National Institute of Science Education and Research, HBNI,
Jatni-752050, India  Poonam Kumari  Department of Condensed Matter Physics and
Material Science, S. N. Bose National Center for Basic Sciences,Kolkata-700106
Priya Mahadevan<EMAIL_ADDRESS>Department of Condensed Matter
Physics and Material Science, S. N. Bose National Center for Basic
Sciences,Kolkata-700106

###### Abstract

The puzzling observation of room temperature ferromagnetism in double
perovskites (A2BB′O6), despite having the magnetic lattice of B-ions diluted
by non-magnetic B′-ions, have been examined for Sr2FeReO6. Ab-initio spin
spiral electronic structure calculations along various high symmetry
directions in reciprocal space are used to determine the exchange interactions
entering an extended Heisenberg model, which is then solved classically using
Monte Carlo simulations to determine the ferromagnetic transition temperature
Tc. We find that one must consider onsite Coulomb interactions at the
nonmagnetic Re sites ($U$) in order to obtain a Tc close to the experimental
value. Analysis of the $ab$-$initio$ electronic structure as well as an
appropriate model Hamiltonian trace the origin of enhancement in Tc with $U$
to the enhanced exchange splitting that is introduced at these sites. This in
turn destabilizes the antiferromagnetic exchange channels, thereby enhancing
the Tc. The role of occupancy at the non-magnetic sites is examined by
contrasting with the case of Sr2FeMoO6.

###### pacs:

##  I Introduction

Double perovskite (DP) oxides have the general formula A2BB′O6 with A being
the rare earth or alkaline-earth ion, B and B′ are the transition metal atoms
separated by oxygen atoms in the lattice. In general, a DP oxide consists of
two sub lattices of perovskite ABO3 and AB′O3 units in a way that the B and B′
sites are arranged alternately in a three dimensional array. B (B′) transition
metal atoms reside in an octahedra formed by six oxygen atoms. The finding of
ferromagnetic behavior in the Re based DP oxides A2FeReO6 (where A= Ba, Sr,
Ca) above room temperature (with ferromagnetic transition temperatures 316,
401, 538 K respectively) ch_sfro_ref12  stimulated research in this new
ferromagnetic class of DP oxides. Another similar example of a high Tc
ferromagnetic DP oxide is Sr2FeMoO6 (SFMO) with a transition temperature of
410 K ch_sfro_ref14 ; ch_sfro_ref14a . In Sr2FeReO6 (SFRO), the separation
between the closest magnetic Fe atoms is $\sqrt{2}$ times larger than that
found in a Fe-based perovskite oxide. This large separation of Fe atoms arises
because of the dilution of the magnetic lattice with non-magnetic Re atoms.
The naive expectation was that it should result in weaker magnetic
interactions and hence a lower magnetic ordering temperature. Surprisingly,
the Curie temperature for this system is 401 K. In contrast, considering the
case of the hole doped manganites manganites_rev , where one has smaller
separations between the magnetic atoms, one finds a lower magnetic ordering
temperature. As the ferromagnetism has been conventionally understood within a
double exchange model, the presence of nonmagnetic atoms like Re/Mo present in
a double perovskite crystal requires a different mechanism to explain the high
ferromagnetic ordering temperature in these compounds.

A qualitative model to understand the high magnetic ordering temperature was
proposed by Sarma et al. ch_sfro_ref16  after analyzing the $ab$ $initio$
electronic structure of SFMO. They concluded that the hopping interactions
between first neighbor Fe and Mo atoms in the lattice led to an exchange
splitting being induced on the Mo sites which is opposite in direction to that
at the Fe site. This led to an antiferromagnetic coupling between the Fe and
Mo sites and an effective ferromagnetic coupling between the Fe sites. This
results in the high magnetic ordering temperatures of these materials.

Electronic structure calculations reveal that the majority Fe-up spin $t_{2g}$
and $e_{g}$ states are $\sim$ 3.75 eV and 2.25 eV below the fermi energy
respectively, while the minority spin states have highly hybridized Fe-Re/Mo
character. This idea led to the use a Kondo-like model to describe the
magnetism in these materials. The majority spin Fe states are approximated by
a spin within this model which interacts with the itinerant conduction
electrons through an antiferromagnetic coupling. The itinerant conduction
electron can delocalize by hopping from the $t_{2g}$ levels on the nonmagnetic
transition metal atom to those on Fe. While this model was able to
successfully explain the high magnetic ordering temperature in Sr2FeMoO6 which
had just one electron on the nonmagnetic atom, it failed to explain the high
magnetic ordering temperature when one had more than one electron on the
nonmagnetic atom as in Sr2FeReO6 ch_sfro_ref17 ; ch_sfro_ref18 ; ch_sfro_ref19
; ch_sfro_ref20 .

Brey and coworkers ch_sfro_ref21  carried out a mean field study considering a
$t_{2g}$ only model. The suppression of Tc was recovered qualitatively by
introducing electron-electron interactions at the nonmagnetic site within the
model. This helped in obtaining the experimental trend of increase in Tc with
increase of the number of electrons at the nonmagnetic site. We reexamine the
changes in the electronic structure at the nonmagnetic site, examining how it
evolves with $U$ and how this affects the stability of the ferromagnetic state
with respect to the competing antiferromagnetic state, in addition to the
modifications in the exchange interaction strengths.

In the context, the issue of high temperature magnetism in Sr2FeReO6 has been
studied within ab initio electronic structure as well as a model Hamiltonian
approach. The $ab$ $initio$ band structure calculations were performed and
electron correlation effects on Fe as well as on Re sites were introduced
within the GGA+U technique. Considering a $U$ of 2.5 eV on the Fe site, the
ground state is found to be ferromagnetic. The magnetic exchange interaction
strengths ($J$’s) between the magnetic (Fe) atoms were calculated using an
inverse Fourier transform method. Using a classical Monte Carlo (MC)
simulation technique employing an extended Heisenberg model, the Tc of the
system was estimated. In the first case, when there was no $U$ on the Re site,
we found that the nearest neighbor exchange interaction $J_{1}$ is
ferromagnetic while the second neighbor exchange interaction $J_{2}$ is
antiferromagnetic contrary to what is believed in literature ch_sfro_ref22 .
$J_{3}$ and $J_{4}$ are negligibly small when compared to $J_{1}$, $J_{2}$.
With these exchange interaction strengths, the estimated Tc within the Monte
Carlo simulation study was much lower ($\sim$ 285 K) than the experimental
value of 401 K. Taking a clue from the work by Brey et al. ch_sfro_ref21 , we
then considered a $U$ of 1 eV on the Re-$d$ states and again calculated the
exchange interaction strengths. We found that the $J_{1}$ got enhanced by a
factor of $\sim$ 2\. Most surprisingly, $J_{2}$ changes its sign from being
antiferromagnetic in the no $U$ case to ferromagnetic in the present case when
a $U$ of 1 eV is considered at the Re site. Using a $U$(Fe) = 2.5 eV and
$U$(Re) = 1.0 eV, we find the magnetic ordering temperature to be 395 K which
is reasonably close to the experimentally observed value. Within $ab$ $initio$
electronic structure calculations, the ferromagnetic ground state is found to
be half-metallic. The Re moments found in this case are induced in nature and
are antiparallel to Fe moments. With an application of $U$ on the Re site, we
found no substantial effect in the density of states for the ferromagnetic
configuration which could explain the gain in stability in this case. However,
examining the density of states in a competing antiferromagnetic configuration
we found that the $U$(Re) played the role on inducing an exchange splitting of
the Re-$d$ states. Thus by destabilizing the antiferromagnetic state, one
found that one could increase the ferromagnetic stability. In order to test
this hypothesis we set up a multiband Hubbard Hamiltonian to describe the
system. Introducing an exchange splitting ($J_{h}$) on Re was found to enhance
the stability of the ferromagnetic state.

##  II Methodology

The electronic and magnetic properties of SFRO have been calculated using a
plane wave pseudopotential implementation of density functional theory using
PAW potentialsch_sfro_ref26 ; ch_sfro_ref27  implemented in the Vienna ab-
initio simulation package (VASP). The generalised gradient approximation (GGA)
ch_sfro_ref28 ; ch_sfro_ref29  has been used for the exchange correlation
functional. The effects of electron-electron correlation on the transitional
metal $d$ electrons were considered within the GGA+U scheme using the Dudarev
et al. formalism ch_sfro_ref30  for the plus U implementation. In our
calculations we have used a $U$ of 2.5 eV on the Fe-$d$ states while a $U$ of
0, 1 and 2 eV was considered on the Re-$d$ states. A k-mesh of 8 $\times$ 8
$\times$ 8 was used for the integrations in momentum space, while a cutoff
energy of 500 eV was used for the plane wave included in the basis for the
calculations. The inter atomic exchange interactions strengths J’s were
determined using the frozen magnon approachch_sfro_ref31 ; ch_sfro_ref32 ;
ch_sfro_ref33 ; ch_sfro_ref34 . In this approach the total energy
corresponding to each spin spiral state associated with a spin spiral vector
$\vec{q}$ was calculated within a noncollinear spin implementation in
VASPch_sfro_ref35 . Each spin at the Fe site is treated classically and
treated as a vector where $\vec{s_{i}}$ is a unit vector along the direction
of the spin at the site labelled by $i$. Hence,

$\vec{s_{i}}=({\sin\theta\cos(\vec{q}\cdot\vec{R_{i}}),\sin\theta\sin(\vec{q}\cdot\vec{R_{i}}),\cos\theta})$
(1)

where $\vec{R_{i}}$ is the position of the atom and $\theta$ is the polar
angle made by the spin with the z-axis. Assuming the spins at the Fe sites to
be localized, we mapped the magnetic energies calculated for various $\vec{q}$
onto an extended Heisenberg model as,

$H=-\frac{1}{2}\sum_{ij}J_{ij}\vec{s_{i}}\cdot\vec{s_{j}}$ (2)

In the above equation, $s_{i}$, $s_{j}$ are unit vectors in the direction of
the spins at sites i and j while $J_{ij}$ are the exchange interaction
strengths between the Fe spins at site $i$ and $j$. The magnitude of spins are
absorbed into the $J$. It is clear that in this model, positive and negative
$J$ value imply ferromagnetic and antiferromagnetic exchange interactions
respectively. Hence, corresponding to each $\vec{q}$ value, one has the total
energy E($\vec{q}$) in terms of the J($\vec{q}$) which is related to the
exchange interactions via Fourier transformation as,

$J(\vec{q})=\sum_{i\neq 0}J_{0i}e^{i\vec{q}\cdot\vec{R_{0i}}}$ (3)

For simplicity, in our calculations, we restricted the spins to lie in the xy
plane After performing an inverse Fourier transform, we obtain the real space
exchange parameters as,

$J_{0i}=\frac{1}{N}\sum_{\vec{q}}J(\vec{q})e^{-i\vec{q}\cdot\vec{R_{0i}}}$ (4)

Thus the calculated spin wave spectra along various various high symmetry
directions with varying $\vec{q}$ has been determined and the spectrum was
fitted using a least squared error minimization procedure to estimate the
exchange parameters with the best fit. The Boltzmann weight factor
$e^{-\beta\Delta E}$ was considered while estimating the exchange parameters
to take care of the fact that statistical weight corresponding to large angles
between the spins are low in the frozen magnon states. Here $\beta=1/k_{B}T$
where $T$ is the ferromagnetic transition temperature of 401 K for SFRO and
$\Delta E=E(\vec{q})-E(0)$ thesis_ashis . This analysis was done for both the
case when there was no $U$ on the Re sites and also with a $U$(Re) = 1 eV
while the $U$ on Fe sites was kept fixed at 2.5 eV.

With the help of these exchange interactions, we have calculated the
ferromagnetic transition temperature (Tc) using a classical Monte Carlo
(MC)ch_sfro_ref36  method. In this technique, we first set up the magnetic
lattice. In SFRO the Fe atoms form a face centered lattice. A system size
consisting of 16 $\times$ 16 $\times$ 16 spins is considered. The extracted
$J$’s from the procedure discussed above are used as inputs for an extended
Heisenberg model. During a single MC run, random orientations ($\vec{s_{i}}$)
for a spin at $i$-th site were considered and the energy of the system was
calculated. Then a new random orientation of spin ($\vec{s^{\prime}_{i}}$) was
created and the system energy was again calculated. If the energy associated
with this new random orientation ($\vec{s^{\prime}_{i}}$) was lower than that
of ($\vec{s_{i}}$), then ($\vec{s^{\prime}_{i}}$) was accepted. Else,
Metropolis algorithmch_sfro_ref37  was employed to decide whether to accept or
reject ($\vec{e^{\prime}_{i}}$). This procedure was performed over all the
lattice sites $i$ is one Monte Carlo step. The magnetization M(T) per spin was
calculated only after bringing the system into thermal equilibrium. Total 5
$\times$ 105 Monte Carlo steps were used in these calculations and results
were found to be well converged with respect to the system size as well as the
MC steps.

The multiband Hubbard Hamiltonian considered for the microscopic understanding
of the system is given as,

$\displaystyle H$
$\displaystyle=\sum_{i,l,\sigma}\epsilon^{f}_{d}d^{\dagger}_{il\sigma}d_{il\sigma}+\sum_{i,l,\sigma}\epsilon^{m}_{d}d^{\dagger}_{il\sigma}d_{il\sigma}+\sum_{i,l,\sigma}\epsilon_{p}p^{\dagger}_{il\sigma}p_{il\sigma}+\sum_{i,l,\sigma}\epsilon_{s}s^{\dagger}_{il\sigma}s_{il\sigma}$
$\displaystyle-\sum_{i,j,l_{1},l_{2},\sigma}\left(t^{l_{1}l_{2}}_{i,j,pd_{f}}d^{\dagger}_{il_{1}\sigma}p_{jl_{2}\sigma}+{\bf
H.c.}\right)-\sum_{i,j,l_{1},l_{2},\sigma}\left(t^{l_{1}l_{2}}_{i,j,pd_{m}}d^{\dagger}_{il_{1}\sigma}p_{jl_{2}\sigma}+{\bf
H.c.}\right)$
$\displaystyle-\sum_{i,j,l_{1},l_{2},\sigma}\left(t^{l_{1}l_{2}}_{i,j,pp}p^{\dagger}_{il_{1}\sigma}p_{jl_{2}\sigma}+{\bf
H.c.}\right)-\sum_{i,j,l_{1},l_{2},\sigma}\left(t^{l_{1}l_{2}}_{i,j,sd_{f}}d^{\dagger}_{il_{1}\sigma}s_{jl_{2}\sigma}+{\bf
H.c.}\right)$
$\displaystyle-\sum_{i,j,l_{1},l_{2},\sigma}\left(t^{l_{1}l_{2}}_{i,j,sd_{m}}d^{\dagger}_{il_{1}\sigma}s_{jl_{2}\sigma}+{\bf
H.c.}\right)$
$\displaystyle+\sum_{\alpha\beta\gamma\delta,\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}U^{\alpha\beta\gamma\delta}_{dd,f}d^{\dagger}_{\alpha\sigma_{1}}d^{\dagger}_{\beta\sigma_{2}}d_{\gamma\sigma_{3}}d_{\delta\sigma_{4}}$
$\displaystyle+\sum_{\alpha\beta\gamma\delta,\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}U^{\alpha\beta\gamma\delta}_{dd,m}d^{\dagger}_{\alpha\sigma_{1}}d^{\dagger}_{\beta\sigma_{2}}d_{\gamma\sigma_{3}}d_{\delta\sigma_{4}}$

where $d^{\dagger}_{il\sigma}\left(d_{il\sigma}\right)$ creates (annihilates)
an electron with spin $\sigma$ in the $l$th $d$ orbital on transition metal
site in the $i$th unit cell while
$p^{\dagger}_{il\sigma}\left(p_{il\sigma}\right)$/$s^{\dagger}_{il\sigma}\left(s_{il\sigma}\right)$
creates (annihilates) an electron with spin $\sigma$ in the $l$th $p$/$s$
orbital on oxygen atom in the $i$th unit cell. Here the index $f$ and $m$
represent the Fe and Re sites respectively. The hopping interaction strengths
($t$) were kept fixed at the values obtained from fitting the nonmagnetic $ab$
$initio$ band structure of Sr2FeReO6 within a tight binding model having Fe
$d$, Re $d$ and O $s$ and $p$ states in the basis. These hopping interactions
were parametrized in terms of Slater-Koster parametersch_sfro_ref38
dp$\sigma$, dp$\pi$, pp$\sigma$, pp$\pi$ and ds$\sigma$. The onsite energies
were however varied, though the same values were taken for both spin channels.
The Coulomb matrix elements were parameterized in terms of the Slater-Condon
integralsch_sfro_ref25  $F^{0}$, $F^{2}$ and $F^{4}$ exploiting its
rotationally invariant nature and there values were determined from atomic
Hartree-Fock calculations. As the value of $F^{0}$ is decreased significantly
from its atomic value due to screening, it was varied and used to determine
the multiplet averaged Coulomb interaction $U$. As the Hund’s exchange
interaction strength is related to the Slater Condon integrals by the relation
$J_{h}=\frac{F^{2}+F^{4}}{14}$, a reduction factor was used to get the desired
exchange splitting observed in the $ab-initio$ spin polarized calculations.
For Fe $d$ orbitals,a value of $J_{h}$ of 0.8 eV was considered while at the
Re site it is varied from 0.24 eV to 0.4 eV.

The four fermion Coulomb interaction terms in the above Hamiltonian were
reduced to quadratic terms assuming a mean field approximation ch_sfro_ref39 ;
ch_sfro_ref40 ; prb4d5d . The expectation values of the order parameters
entering the mean field decoupled Hamiltonian were determined self-
consistently till an energy convergence of 10-5 eV was achieved. Total energy
and density of states were calculated using the tetrahedron method of
integration for a 10 $\times$ 10 $\times$ 10 k-mesh.

##  III Results and Discussion

An ordered double perovskite structure consists of two different transition
metal atoms B and B′ arranged alternately, each with a local environment
involving BO6 and B′O6 octahedra. The crystal structure of Sr2FeReO6 is simple
cubic with space group Fm-3m sfro_exp . We have optimised both the unit cell
parameters as well as the internal positions within ab-initio electronic
structure calculations with $U$(Fe) = 2.5 eV. The theoretically optimized
structural information are given in Table 1.

Structural details | Experimental | Experimental | Theoretical
---|---|---|---
Space group | Fm-3m | I 4/m | Fm-3m
a = b (Å) | 7.888 (= c) | 5.561 | 7.890 (= c)
c (Å) | – | 7.901 | –
$\alpha$ = $\beta$ = $\gamma$ (°) | 90 | 90 | 90
Bond lengths (Å) |  |  |
Re-O | 1.70 | 1.959/1.962 | 1.94
Fe-O | 2.24 | 1.985/1.991 | 2.01
Fe-O-Mo Bond angles (in degrees) |  |  |
In-plane | 180 | 170 | 180
Out-of-plane | 180 | 180 | 180
Cell volume (Å3) | 490.8 | 488.8 | 491.17
Table 1: Experimental and optimized/theoretical structural details for
Sr2FeReO6. In second and third column, experimental details given for cubic
sfro_exp  and tetragonal sfro_tetra  structures respectively. Fe-O and Re-O
bond lengths are quite similar for the experimental tetragonal phase and our
optimized crystal structure of Sr2FeReO6. phase

The experimental unit cell parameters are also given for comparison. As one
can see, in SFRO all the Fe-O-Re bond lengths are 180° while all the Fe-O or
Re-O bondlengths are identical within the respective octahedra as is evident
from Table 1. Further there was no substantial change in the experimentally
observed lattice constants found after structural optimization. However, there
is substantial change in the Fe-O and Re-O bond lengths after the
optimization. We use this optimized structure in the rest of our calculations.
Obtained Fe-O and Re-O bond lengths in for our optimized structure is quite
similar with the ones found in tetragonal phase of Sr2FeReO6. This similarity
among the bond lengths of the tetragonal crystal structure along with the fact
that the tetragonal distortion is quite small in magnitude ($\sim$ 0.011 Å)
justify use of optimized crystal structure in our calculations. To understand
the high magnetic ordering temperature of this system, the magnon dispersions
were calculated along various symmetry directions in the Brillouin zone of
SFRO (Fig. 1). Each point in Fig. 1 corresponds to a different magnetic
configuration identified by $\vec{q}$ and has been calculated for a $U$ of 2.5
eV on the Fe $d$ states. Clearly, $\vec{q}$=0, ($\Gamma$) which represents the
ferromagnetic state has the lowest energy.

Figure 1: Magnon dispersion plot for Sr2FeReO6 for $U$(Fe) = 2.5 eV and
$U$(Re) = 0 and 1 eV. The G-point energy is the ferromagnetic ground state
energy which gets higher stability with application of $U$ = 1 eV on Re site.

This is in agreement with the experimental findings. The value of magnetic
moment at Fe and Re sites are 3.9 $\mu_{B}$ and -0.95 $\mu_{B}$ respectively.
The induced moment on the Re site is opposite to the Fe moments and originates
from Fe-Re hopping as described in Refch_sfro_ref16 . This is evident from the
density of states calculated for a perfectly antiferromagnetic configuration.
In that case one had no exchange splitting of the Re-$d$ states, which one
would otherwise find if there was a component emerging from an intrinsic
exchange splitting.

Figure 2:  Fe-Fe interatomic exchange interaction parameters ($J$) of
Sr2FeReO6 as a function Fe-Fe nearest neighbor distance (in Å).

The real space exchange interaction strengths ($J$’s) are then estimated from
an inverse Fourier transformation. A plot of the different $J$’s as a function
of distance is given in Fig. 2. We found that in this case when there is no
$U$ on the Re-$d$ states, $J_{1}$ which is the most dominant one was found to
be 10.9 meV and is ferromagnetic while $J_{2}$ is -5.1 meV and is
antiferromagnetic in nature. All other J’s are negligibly small. A schematic
showing the first four exchange interactions and the pathways they represent
is shown in Fig 3. One should mention here that the separation between the Fe
atoms connected via $J_{1}$, $J_{2}$, $J_{3}$ and $J_{4}$ are $\sim$ 5.58 Å,
7.89 Å, 9.66 Å and 11.16 Å respectively. Using these values for the exchange
interaction strengths in a Monte Carlo calculation, one finds a magnetic
transition temperature of 285 K. As shown by Brey et al., application of $U$
on the nonmagnetic Mo site will penalize double occupancy. It consequently
leads to the ferromagnetic state being favored. Taking a clue from this
hypothesis, we applied a $U$ of 1 eV on the Re-$d$ states while a $U$ on
Fe-$d$ states was kept fixed at the value of 2.5 eV and recalculated the
magnon spectra shown in Fig 1 (red dashed lines). A plot of estimated values
of J’s in this case is shown in Fig 2.

Figure 3: A schematic to show the first four Fe-Fe interatomic exchange
interaction strengths ($J$’s) in the simple cubic lattice of Sr2FeReO6.

A clear enhancement in the stability of the ferromagnetic configuration can be
found in this case when compared to the previous case of no $U$ on Re. The
effect of $U$ on the Re site is also manifested in the exchange parameters. An
analysis of the estimated $J$’s in this case yields mainly two important
points. First is the enhancement in the magnitude of the most dominant nearest
neighbor ferromagnetic exchange interaction strength $J_{1}$. Its value $\sim$
21 meV is nearly doubled with a $U$(Re) = 1 eV from its value in the absence
of $U$ on Re. Most interestingly, $J_{2}$ which was antiferromagnetic in the
previous case is now found to be ferromagnetic in the present case. The
estimated Tc with this new set of exchange interaction strengths is 395 K
which is in reasonable agreement with the experimental value of 401 K. The
magnetic moment on the Re site is enhanced to -1.04 $\mu_{B}$ in this case
while it is the same for the Fe atoms when there was no $U$ on the Re-d
states. A plot of the normalized magnetic moment variation as a function of
temperature in both the cases when there is no $U$ on Re and when a $U$ of 1
eV is considered on the Re-$d$ states is shown in Fig 4. There is indeed a
clear enhancement of the Tc.

Figure 4:  Variation of normalized magnetization as a function of temperature
calculated for $U$(Fe) = 2.5 eV and $U$(Re) = 0 and 1 eV using Monte Carlo
simulations for Sr2FeReO6.

To explain the role of a $U$ on the Re $d$ states, we went to examine the Fe
and Re atom projected density of states in both the cases with and without a
$U$ on Re-$d$ states for the ferromagnetic ground states, a plot of which is
shown in left panel of Fig. 5. We also include the case where a $U$ of 2 eV
was considered on the Re $d$ states. The large exchange splitting at the Fe
site and an induced small exchange splitting at Re site is evident from the
plot.

Figure 5: Atom projected Re-$d$ density of states (DOS) plots from $ab$
$initio$ calculations for Sr2FeReO6. The plot on the left is for a
ferromagnetic configuration of Fe spins while the plot in the right is from
the closest competing antiferromagnetic state of Fe spins. The zero of energy
is the Fermi energy.

As we can see in the up spin channel, Fe states are fully occupied verifying
the 3$d^{5}$ ionic configuration of Fe while the Re states are unoccupied and
are $\sim$ 0.6 eV above the Fermi level. With the application of $U$ on Re,
the unoccupied Re-up spin states are pushed up away from the Fermi level. A
$U$ of 1 eV on Re pushes the Re $d$ states by $\sim$ 0.25 eV from its position
for the case when there was no $U$. In the down spin channel, there is hardly
any significant change in Fe and Re-$d$ density of states with the variation
of $U$ on the Re site. They are found to be pinned to their positions in the
absence of $U$. Thus, the small shifts in the Re-up spin density of states
above the Fermi energy as function of $U$ on the Re $d$ states alone is
therefore insufficient to describe the substantial gain in the stability of
the ferromagnetic ground state found by examining the magnon dispersion shown
in Fig. 1. For example, the gain is $\sim$ 50 meV at the X and W points while
it is even more than $\sim$ 85 meV at the L-point over the $U$ = 0 values. We
went on to examine the density of states of the closest competing
antiferromagnetic state. The spins on the Re site are antiferromagnetically
coupled to the spins on the Fe sites. We consider antiferromagnetic
configurations in which we have flipped the spin of two out of the six Fe
nearest neighbors of a Re atom and made them ferromagnetic. As this is a
collinear antiferromagnetic configuration with the smallest number of spins
around Re flipped, one would expect this to be the closest lying
antiferromagnetic excited states.

The atom projected density of states for Fe and Re atoms are then calculated
by varying the $U$ on the Re site. The calculated density of states for the
Fe-$d$ states have almost similar features as found for the ferromagnetic
configuration in Fig. 5. Similar to the ferromagnetic case, these states are
unaffected by the application of $U$ on Re atom and hence are not shown here.
It is the Re-$d$ density of states shown in the right panel of Fig. 5, which
have the most dramatic effects in this antiferromagnetic configuration with
application of $U$ on it. When there is no $U$ on the Re states, both spin
channels are found to be contributing at the Fermi level as shown in the right
panel (a) of Fig. 5. The induced exchange splitting at the Re site, which was
present in the ferromagnetic case has reduced substantially. This provides
pathways for both, the ferromagnetically coupled as well as the
antiferromagnetically coupled spins. This explains the case of reduced
magnetic transition temperature. The Re-up spin channel gets substantially
depleted at the Fermi level with the application of $U$ on Re-$d$ states as
shown in the right panel (b) of Fig. 5. The majority weight of the Re-up spin
states is now much above the Fermi level bringing an effective exchange
splitting at the Re site. Because of this the energy gain from delocalization
of antiferromagnetically coupled Fe spins is significantly reduced. The
exchange splitting at the Re site, further enhances the energy gain for
ferromagnetically coupled spins. This is manifested in terms of enhanced
$J_{1}$ and a ferromagnetic $J_{2}$ in this case and consequently a higher Tc.
This picture gets even more clear by looking at Re density of states with a
$U$ of 2 eV at its $d$ states where the effects are more pronounced. A
comparison of the stability of the ferromagnetic states with respect to this
antiferromagnetic state as a function of $U$(Re) is tabulated in Table 2.

AFM-FM (eV)
---
URe (eV) | Sr2FeReO6 | Sr2FeMoO6
0.0 | 0.085 | 0.129
1.0 | 0.138 | 0.143
2.0 | 0.175 | 0.158
Table 2: Change in stability of the ferromagnetic state calculated with
respect to closest competing antiferromagnetic state from $ab$ $initio$
calculations as function of $U$ on the Re/Mo-$d$ states for Sr2FeReO6 and
Sr2FeMoO6 double perovskites.

In order to explicitly understand the role of the exchange splitting at the
nonmagnetic Re site, we further went on to develop a simple model which can
define this system correctly. For this purpose, we start with calculating the
$ab$ $initio$ atom projected density of states within a non-magnetic
calculation for Sr2FeReO6 and is given in Fig. 6. The Fe $d$ states are found
in the energy window from -2 eV to 1 eV. The octahedral field of the six
oxygens which are the nearest neighbors of Fe lift the degeneracy of the $d$
orbitals and one has the $t_{2g}$ orbitals contributing in the energy window
from -2 to -0.5 eV. The $e_{g}$ states contribute in the energy window from
$\sim$ 0 to 1 eV. We find a similar degeneracy lifting of the $d$ states on
Re, with the $t_{2g}$ states contributing in the energy window from $\sim$ 0
to 0.5 eV. Examining the O $p$ states, one finds that their primary
contribution is below $\sim$ -3 eV and extends to -8 eV below the Fermi level.
There is a small admixture of the O $p$ states in the energy window
contributed primarily by the Fe and Re $d$ states. This arises from $p$-$d$
hybridization and one finds that the contributions are larger for the Fe
$e_{g}$ states as well as the Re $t_{2g}$ states which have the larger hopping
matrix elements associated with them. Surprisingly we find the energy window
over which the Fe $t_{2g}$ states contribute comparable with that in which Fe
$e_{g}$ states contribute. This seems surprising and could emerge from the
smaller energy denominator for the delocalization of an electron in the Fe
$t_{2g}$ states via the Re $t_{2g}$ states than the Fe $e_{g}$ states via the
levels with the same symmetry on Re.

Figure 6: Fe-$d$, Re-$d$ and oxygen-$p$ projected partial density of states
for non-magnetic Sr2FeReO6 from $ab$ $initio$ calculations. The zero of energy
is the Fermi energy. Figure 7:  Non-magnetic band dispersions from (a) $ab$
$initio$ , (b) tight binding calculations with Fe-$d$, Re-$d$ and O-$s$ and
$p$ basis states along various symmetry directions in Brillouin zone of
Sr2FeReO6. The zero of energy is the Fermi energy.

The ab initio band dispersions along various symmetry directions for
nonmagnetic Sr2FeReO6 are shown in Fig. 7 in the energy window -4 to 2 eV. As
Fe $d$, Re $d$ and O $p$ states were found to contribute in this energy
window, we included these states in the basis of tight binding model
considered. We also included the $s$ states on the oxygen as discussed in the
methodology section. The parameters entering the tight binding Hamiltonian
should reproduce the band structure from the $ab$ $initio$ calculations. A
least square error minimization procedure was used. The best fit parameters
are given in Table 3.

Onsite energies |  |  | Fe-O |  | Re-O |  | O-O
---|---|---|---|---|---|---|---
Fe-$d$ | Re-$d$ | O-$p$ |  | $dp\sigma$ | $dp\pi$ | $ds\sigma$ |  | $dp\sigma$ | $dp\pi$ | $ds\sigma$ |  | $pp\sigma$ | $pp\pi$ |
$t_{2g}$ | $e_{g}$ | $t_{2g}$ |  |  |  |  |  |  |  |  |  |  |  |  |
2.46 | 3.34 | 3.90 | 0.68 |  | 1.90 | -1.495 | -1.44 |  | 3.05 | -1.28 | -0.76 |  | 0.75 | -0.164 |
Table 3: Onsite energies and various hopping interaction strengths (in eV)
parameterized in terms of Slater-Koster parameters obtained by fitting $ab$
$initio$ band structure of non-magnetic Sr2FeReO6 employing a tight binding
model with Fe, Re -$d$ and O $s$ and $p$ states in the basis.

The $dp\pi$ interactions between the Fe $d$ and O $p$ states are significantly
enhanced compared to a ratio of 0.5 for $dp\sigma/dp\pi$ interactions that one
usually expects explaining partially the similar bandwidths that we find for
Fe $t_{2g}$ and Fe $e_{g}$ states. $ds\sigma$ interactions control the crystal
field splitting between the $t_{2g}$ and $e_{g}$ orbitals and are found to be
considerable for the Fe $d$ states also, compared to the Re $d$ states as the
Re $e_{g}$ states were not included in the fitting. Examining the onsite
energies, one finds that one has the Fe $t_{2g}$ states, Fe $e_{g}$ states and
then the Re $t_{2g}$ states as expected.

Onsite energies |  | Fe-Fe |  | Fe-Re |  | Re-Re
---|---|---|---|---|---|---
Fe-$d$ | Re-$d$ |  | $dd\sigma$ | $dd\pi$ | $dd\delta$ |  | $dd\sigma$ | $dd\pi$ | $dd\delta$ |  | $dd\sigma$ | $dd\pi$ | $dd\delta$
$t_{2g}$ | $e_{g}$ | $t_{2g}$ |  |  |  |  |  |  |  |  |  |  |  |
4.45 | 6.39 | 5.10 |  | -0.05 | 0.05 | -0.004 |  | -0.866 | 0.298 | 0.0 |  | -0.199 | 0.008 | -0.02
Table 4: Onsite energies and various hopping interaction strengths (in eV)
parameterized in terms of Slater-Koster parameters obtained by fitting $ab$
$initio$ band structure of non-magnetic Sr2FeReO6 employing a tight binding
model with only Fe and Re -$d$ states in the basis.

The charge transfer energy between the O $p$ states and Fe $t_{2g}$ states is
small but positive. In order to examine the role of explicitly including the O
$p$ states in the basis, especially when the charge transfer energy is so
small, a mapping was carried out onto a $d$-only model. The parameters found
to give the best fit of the $ab$ $initio$ band structure shown in Fig 8(a) are
given in Table 4. The tight binding band structure calculated for the best fit
parameters are given in Fig. 8(b).

Figure 8: Non-magnetic band structure from (a) $ab$ $initio$ , (b) tight
binding calculations with only Fe and Re $d$ basis states, plotted along
various symmetry directions in Brillouin zone of Sr2FeReO6. The zero of energy
is the Fermi energy.

As is evident, we include explicitly the Fe $d$ and Re $d$ states while
constructing the tight binding Hamiltonian with $d$-$d$ interactions.
Interestingly we find that the Re $t_{2g}$ levels lie in between the Fe
$t_{2g}$ and Fe $e_{g}$ levels. This is because in the effective $d$-only
model, we are approximating the positions of the $d$ levels to their positions
after their interaction with oxygen. The matrix elements for Fe $e_{g}$ and O
$p$ interactions are larger, resulting in larger shifts compared to the
interaction between Re $t_{2g}$ and O $p$ states. This leads to negative
energy difference between the Fe $e_{g}$ and the Re $t_{2g}$ states. In order
to verify if the placing of the Re $t_{2g}$ levels below the Fe $e_{g}$ levels
is not an artefact of the fitting procedure, a similar analysis was carried
out for a $d$ only mapping using maximally localized wannier
functionsch_sfro_ref44  for the radial parts of the basis functions. A
comparison of the $ab$ $initio$ band structure and the fitted band structure
are given in Fig. 9.

Figure 9:  Non-magnetic $ab$ $initio$ band structure fitted with a Wannier
functions based tight binding model considering Fe and Re -$d$ orbitals in the
basis. The zero of energy is the Fermi energy.

Here we find that the fit is reasonable which gives us confidence in examining
the onsite energies. We find that the difference between the Re $t_{2g}$ and
the Fe $e_{g}$ levels is negative, $\sim$ -0.7 eV in agreement with what we
had earlier. In order to determine the appropriate model, a multiband Hubbard
model is considered with the same value of $U$ and the electronic structure
for a given set of parameters is calculated using a mean field decoupling
scheme for the four Fermion terms. The hopping interaction strengths are kept
fixed at the values obtained by fitting the nonmagnetic band structure.

We consider the charge fluctuations at the Fe and Re sites as,

$\displaystyle\Delta_{Fe\leftrightarrow
O}=E_{f}(d^{6}\bar{L})-E_{i}(d^{5}L^{0})=\epsilon_{d}-\epsilon_{p}+5U^{Fe}$
(5) $\displaystyle\Delta_{Re\leftrightarrow
O}=E_{f}(d^{3}\bar{L})-E_{i}(d^{2}L^{0})=\epsilon_{d}-\epsilon_{p}+2U^{Re}$
(6)

Here $\epsilon_{d}$ and $\epsilon_{p}$ are the onsite energies of Fe/Re-$d$
and O-$p$ orbitals respectively. Hence $\Delta$ was kept fixed with
$U_{Fe/Re}$ by changing the onsite energies in accordance to the above
equations. The $J_{h}$ are varied till a good description of the $ab$ $initio$
density of states is obtained. A comparison is shown in Fig 10.

Figure 10:  Atom projected density of states plots from $ab$ $initio$
calculations (upper two) and model Hamiltonian calculation (lower two) with Fe
and Re $d$ states in the basis of Sr2FeReO6. The zero of energy is the Fermi
energy.

A $J_{h}$ of 0.8 eV was considered on Fe and this gave a good description of
the $ab$ $initio$ density of states. Additionally the charge transfer energy
between both Fe and Re with oxygen were found to be 0.85 and 3.5 eV, both
positive. Additionally a $U$ of 0.26 eV and a small $J$ of 0.24 eV was used on
the Re states to get the half-metallic character of the density of states. As
discussed previously with the $ab$ $initio$ results, a $U$ on the Re-$d$
states introduces an effective exchange splitting on the same, resulting in a
gain for the ferromagnetic ground state. We probe this fact independently
using the model Hamiltonian that we have developed. We calculate the stability
of the ferromagnetic ground state with respect to the previously considered
antiferromagnetic state with variation of exchange splitting ($J_{h}$) at the
Re site.

The calculations were done with a fixed $U$ of 0.4 eV of the Re-$d$ states and
keeping all the other parameters the same as their previous values. The
obtained results are tabulated in Table 5.

Jh (eV) | AFM-FM (eV)
---|---
0.24 | 0.134
0.35 | 0.173
0.50 | 0.186
Table 5: Within a multiband Hubbard model, obtained change in stability of the
ferromagnetic state with respect to closest competing antiferromagnetic state
as function of exchange splitting at the Re site for Sr2FeReO6.

As we can see that with an increase in $J_{h}$ from 0.24 eV to 0.35 eV at the
Re site, the stability of the ferromagnetic state is indeed enhanced by $\sim$
40 meV. This supports our previously given explanation for increase in the
stability of ferromagnetic state with application of $U$ on the Re site as it
causes an effective exchange splitting at the Re site. This suggests a route
to enhance the Tc of these double perovskite systems by using magnetic atoms
at the B′ site.

Figure 11: Atom projected Fe, Mo-$d$ density of states plots from $ab$
$initio$ calculations for Sr2FeMoO6. The plot in the left is obtained for a
ferromagnetic configuration of Fe spins while the plot in the right is from
the closest competing antiferromagnetic state. The zero of energy is the Fermi
energy.

Finally, a comparison with a similar system Sr2FeMoO6 (SFMO) has been drawn to
explain the role two of electrons present in the Re-$d$ orbitals (Re5+) in the
case of SFRO. Change in the ferromagnetic stability with respect to the same
antiferromagnetic configuration is calculated for SFMO and the results are
tabulated in Table 2. As one can see that in the case of SFMO, the stability
of ferromagnetic state increases when a $U$ of 1 eV is applied at the Mo site,
but the gain is much less than that of SFRO. While in the case of SFRO this
gain was $\sim$ 50 meV, here it is only 14 meV. Further increase of $U$ at the
Re site from 1 eV to 2 eV causes the stability to decrease in the case of SFRO
while the stability is nearly constant in the case of SFMO. This weaker effect
of $U$ on the Mo site in bringing an effective exchange splitting on it can
easily be seen by analyzing the density of states for Fe, Mo-$d$ orbitals. The
features in the density of states shown in Fig.11 for the ferromagnetic
configuration (left panel) is almost similar to that of the SFRO with the only
difference that the weight of Mo-$d$ orbitals in SFMO is less at the Fermi
level when compared to Re-$d$ states in SFRO. This can be understood as the Re
atom has two electrons in its $d$ orbitals in the ionic configuration of 5+
while Mo has only one electron in its 4$d$ orbital. In the antiferromagnetic
configuration, $U$ of 1 eV on the Mo-$d$ has almost no effect on the its
density of effects while a $U$ of 2 eV brings only a small change at the Fermi
level explaining the less gain in stability of the ferromagnetic state with
the variation of $U$.

##  IV Conclusions

In this study, using $ab$ $initio$ electronic structure calculations we have
explained the role of electron-electron correlations at the nonmagnetic Re
site in Sr2FeReO6 in the context of its high ferromagnetic ordering
temperature. We show that the role of $U$ on Re site is to introduce an
effective exchange splitting reducing the gain from the antiferromagnetic
pathways, providing higher stability to the ferromagnetic ground state. A
simple multiband Hubbard like model is developed to independently verify this
finding. Within this model, indeed, we find an increase in the stability of
the ferromagnetic ground state with increase of exchange splitting at Re site.
Role of two electrons of Re-$d$ orbitals result in substantial gain of the
ferromagnetic stability and this is also contrasted with a similar compound,
Sr2FeMoO6, which has single $d$ electron at the Mo site.

## References

  * (1) J. Longo and R. Ward, J. Am. Chem. Soc. 83, 2816 (1961).
  * (2) K.-I. Kobayashi, T. Kimura, H. Sawada, K. Terakura and Y. Tokura, Nature 395, 677 (1998).
  * (3) K.-I. Kobayashi, T. Kimura, Y. Tomioka, H. Sawada, K. Terakura and Y. Tokura, Phys. Rev B 59, 11159 (1999).
  * (4) J. M. D. Coey, M. Viret, and S. Von Molnar, Adv. Phys. 48, 167 (1999).
  * (5) D. D. Sarma, P. Mahadevan, T. Saha-Dasgupta, S. Ray and A. Kumar, Phys. Rev. Lett. 85, 2549 (2000).
  * (6) J. Kondo, Prog. Theor, Phys. 32, 37 (1964).
  * (7) K. Phillips, A. Chattopadhyay and A. J. Millis, Phys. Rev. B 67, 125119 (2003).
  * (8) J. L. Alonso, and L. A. Fernández, and F. Guinea, and F. Lesmes and V. Martín-Mayor’ Phys. Rev. B 67 214423 (2003).
  * (9) G. Jackeli, Phys. Rev. B 68, 092401 (2003).
  * (10) L. Brey, M. J. Calderon, S. Das Sarma and F. Guinea, Phys. Rev. B 74, 094429 (2006).
  * (11) O. Erten, O. N. Meetei, A. Mukherjee, M. Randeria, N. Trivedi and P. Woodward, Phys. Rev. Lett. 107, 257201 (2011).
  * (12) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * (13) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * (14) G. Kresse and J. J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * (15) J. P. Perdew, K. Burke and Y. Wang, Phys. Rev. B 54, 16533 (1996).
  * (16) S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys and A. P. Sutton, Phys. Rev. B 57, 1505 (1998).
  * (17) L. M. Sandratskii, physica status solidi (b) 136 (1), 167 (1986).
  * (18) S. V. Halilov, H. Eschrig, A. Y. Perlov and P. M. Oppeneer, Phys. Rev. B 58, 293 (1998).
  * (19) O. Grotheer, C. Ederer, M. Fähnle, Phys. Rev. B 63, 100401 (2001).
  * (20) Ležaić, Marjana and Mavropoulos, Phivos and Bihlmayer, Gustav and Blügel, Phys. Rev. B 88, 134403 (2013).
  * (21) M. Marsman and J. Hafner, Phys. Rev. B 66, 224409 (2002).
  * (22) P. Mahadevan, N. Shanthi and D. D. Sarma, Phys. Rev. B 54, 11199 (1996).
  * (23) J. C. Slater and G. F. Koster, Phys. Rev. 94, 1498 (1954).
  * (24) J. C. Slater, Phys. Rev. 34, 1293 (1929); E. U. Condon, Phys. Rev. 36, 1121 (1930); see also, J. S. Griffith, The Theory of Transition-Metal Ions, Cambridge University Press, Cambridge, (1961).
  * (25) D. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics, Cambridge University Press, Cambridge, England, (2000).
  * (26) N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller and E. Teller, J. Chem. Phys. 21, 1087 (1953).
  * (27) J. C. Slater and G. F. Koster, Phys. Rev. 94, 1498 (1954).
  * (28) S. Nimkar, D. D. Sarma, H. R. Krishnamurthy, and S. Ramasesha, Phys. Rev. B 48, 7355 (1993).
  * (29) P. Mahadevan, K. Sheshadri, D. D. Sarma, H. R. Krishnamurthy and R. Pandit, Phys. Rev. B 55, 9203 (1997).
  * (30) A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt and N. Marzari, Comp. Phys. Comm. 178 685 (2008).
  * (31) Ashis Kumar Nandy, Ph.D. Thesis, University of Calcutta (2013).
  * (32) S. Middey, Ashis Kumar Nandy, S. K. Pandey, Priya Mahadevan and D. D. Sarma, Phys. Rev. B 86, 104406 (2012).
  * (33) A. W. Sleight, J. Longo and R. Ward, Inorg. Chem. 1 245 (1962).
  * (34) H. Jalili, N. F. Heinig and K. T. Leung, Phys. Rev. B 79 174427 (2009).
  * (35) H. Kato, T. Okuda, Y. Okimoto, Y. Tomioka, K. Oikawa, T. Kamiyama and Y. Tokura, Phys. Rev. B 69 184412 (2004)