text:
# Simplicial edge representation of protein structures and alpha contact potential with confidence measure

## 1 Introduction

Potential function plays important roles for a variety of computational studies of proteins, including prediction of protein structures, characterization of ensemble thermodynamic properties of proteins, and design of novel proteins. For example, an essential requirement for the prediction of three-dimensional structure of protein from primary sequence is a potential function which can select the native conformation from an ensemble of alternative conformations. Potential function is also often used to guide the sampling of protein conformations . A variety of potential functions have been developed for these important tasks, including physical model based potentials , empirical statistical potentials , and empirical potentials obtained from optimization .

The effectiveness of potential function depends on another critically important factor, i.e., the representation of protein structures. Within this framework, we explore a new type of pairwise contact potentials using a novel contact definition. We introduce a contact definition that reflects protein geometry more accurately. These contacts are based on the computation of the alpha shape, or the dual simplicial complex description of protein structures . Here contact occurs if atoms from non-bonded residues share a Voronoi edge, and this edge is at least partially contained in the body of the molecule. In another word, atoms have to be in physical nearest neighbor contact. In this study, we only examine the 1-simplices, or edges in the dual simplicial complex, which represents the pairwise contacts. This description is related to the work from Wodak’s group, where the contact area between atoms are calculated as the area of intersection of the accessible atom ball around each atom and the faces of its weighted Voronoi cell .

We develop a statistical contact propensities based on alpha edge simplices and a combinatorial null model. To account for the uncertainty of estimated propensity parameter, we develop a studentized bootstrap method for estimating their 95% confidence intervals. We also examine how geometric representation of dual simplicial complex influence the effectiveness of pairwise contact potential functions. In addition, we aim to understand how the size of alphabet of amino acid residues affect the effectiveness of empirical pair potentials. This is important for protein design, where any reduction of the alphabet size of residues will result in exponentially more efficient sampling in the sequence space, therefore leading to more successful design strategies .

This work is also motivated by the need to develop potential functions that takes advantage of recent development in computational geometry and computational topology. The alpha shape representation of proteins from computational geometry allow rapid calculations of metric, topological, and combinatorial structures of proteins precisely . These advantages can lead to improvements in voids and binding site detection , in hierarchical representation of protein dynamic shapes at different resolution, and in conformation sampling. Recent theoretical development suggests many intriguing applications in protein studies . These important advances are largely unexploited currently, and we hope this work provides a useful link by developing empirical pairwise contact potentials based on dual simplicial complex representations of proteins.

Our approach also solves a problem that cannot be satisfactorily addressed with current contact pair potentials. In current approaches, pairwise contact interactions are declared if two residues are within a specific cut-off distance. Potentials based on this contact definition have achieved considerable success. Nevertheless, contacts by distance cut-off can potentially include many implausible non-contacting neighbors, which have no significant physical interaction . Whether or not a pair of residues can make physical contact depends not only on the distance between their center positions (such as Cα or Cβ, or geometric centers of side chain), but also on the size and the orientations of side-chains . Furthermore, two atoms close to each other may in fact be shielded from contact by other atoms. As emphasized in , these contact pairs should not contribute to the assessment of pairwise contacts. By occupying the intervening space, other residues can block a pair of residues from interacting with each other. Inclusion of these fictitious contact interactions would be undesirable.

We organized this paper as follows. First, we describe briefly the dual simplicial complex representation of proteins structures. We then discuss the probabilistic models for developing pairwise potentials. A bootstrap resampling procedure is then introduced which provide confidence intervals of estimated pairwise potential parameters. We then present the pairwise contact potential along with experimental results in discriminating native structure against decoy structures using several benchmark data sets. We further examine how pair potentials developed from the dual simplices compare with cut-off contact definitions using side-chain centers. The effects of reduced alphabet size for amino acid residues are then described. We conclude with remarks and discussion.

## 2 Model and Methods

### Alpha Contacts from Dual Simplicial Complex.

Alpha shape has been successfully applied to study a number of problems in proteins, including void measurement, binding site characterization, protein packing, electrostatic calculations, and protein hydrations . Details of alpha shape have been described elsewhere, here we only provide a brief description for completeness.

To illustrate, Figure 1a shows a two-dimensional molecule formed by a collection of disks of uniform size. Each Voronoi cell is defined by its boundaries, shown as broken lines. Every Voronoi edge is a perpendicular bisector of the line between two atom centers. Each Voronoi cell contains one atom, and every point inside a Voronoi cell is closer to this atom than to any other atom. Three connected Voronoi edges meet at a Voronoi vertex. Another geometric construct, the Delaunay triangulation (Fig 1b) is mathematically dual to the Voronoi diagram, and can be explained by the following procedure. For each Voronoi edge, connect the corresponding two atom centers with a line segment, and for each Voronoi vertex, place a triangle spanning the three atom centers of the three Voronoi cells. Completing this for all Voronoi edges and Voronoi vertices gives a collection of line segments and triangles. Together with the vertices representing atom centers, they form the “Delaunay complex”, which is the underlying structure of Delaunay triangulation.

Now we remove all Delaunay edges (or line segments) where corresponding Voronoi edge of the two atoms does not intersect with the molecular (Figure 1c)). When two atoms are spatially very close, the balls representing the two atoms intersect, and these two atoms have non-zero, two-body volume overlap. When three atoms are spatially very close, they intersect and have non-zero, three-body volume overlap. We further remove all Delaunay triangles where the corresponding Voronoi vertex of the three atoms is not contained within the molecule. The subset of the Delaunay complex formed by the remaining triangles, edges and vertices (atom centers) is called the dual simplicial complex, or the alpha complex. We are interested in identifying only contacting atoms that are spatial nearest neighbors. These are precisely atoms with two-body volume overlaps whose Voronoi cells intersect. By following the mathematical dual structure, i.e., the edges in the alpha complex, we can accurately identify all contacting nearest neighboring atom pairs. For convenience, we use a rather arbitrary criterion and declare two residues are in alpha contact if there is at least one edge connecting these two residues.

Using the alpha shape API kindly provided by Prof. Edelsbrunner, a program interface has been implemented to compute contacting atoms based on precomputed Delaunay triangulation and alpha shape. The Delaunay triangulation is computed using the delcx program, and the alpha shapes computed using the mkalf program. Both can be downloaded from the website at NCSA (http://www.ncsa.uiuc.edu). The van der Waals radii of protein atoms are taken from . We follow Singh & Thornton and increment the van der Waals radii by 0.5 Å. This increment is small and comparable with the resolution of the structure. It enables the modeling of imprecisely determined atomic coordinates without introducing many spurious two-body contacts.

### Probabilistic Model for Pairwise Alpha Contact Propensity.

The propensity $`P(i,j)`$ for residue of type $`i`$ interacting with residue of type $`j`$ is modeled as the odds ratio of the observed probability $`q(i,j)`$ and the expected probability $`p(i,j)`$ of a pairwise alpha contact involving both residue $`i`$ and $`j`$, :

$$P(i,j)=\frac{q(i,j)}{p(i,j)}$$
(1)
To compute $`p(i,j)`$ and $`q(i,j)`$, we choose a simple null model, where the observed contacts of different proteins in the entire database are pooled together, and the expected contact numbers are calculated. This is the same as the reference state of composition-independent scale discussed in literature . We have:

$$q(i,j)=\frac{a(i,j)}{_{i^{},j^{}}a(i^{},j^{})}$$

Here, $`a(i,j)`$ is the number count of atomic contacts between residue type $`i`$ and residue type $`j`$, and $`_{i^{},j^{}}a(i^{},j^{})`$ is the total number of all atomic contacts. The observed probability is then compared against the random probability $`p(i,j)`$ that a pair of contacting atoms is picked from a residue of type $`i`$ and a residue of type $`j`$, when chosen randomly and independently . We have:

$$p(i,j)=N_iN_j(\frac{n_in_j}{n(nn_i)}+\frac{n_in_j}{n(nn_j)}),\text{ when}ij$$

and

$$p(i,j)=N_iN_{i1}\frac{n_in_i}{n(nn_i)},\text{ when }i=j$$

where $`N_i`$ is the number of interacting residues of type $`i`$, $`n_i`$ is the number of atoms residue of type $`i`$ has, and $`n`$ is the total number of interacting atoms.

The alpha contact potential between residue $`i`$ and residue $`j`$ is obtained from the propensity value $`P(i,j)`$ as $`\mathrm{ln}P(i,j)`$, and the overall energy of a protein is calculated as:

$$ϵ=\underset{i,j}{}\mathrm{ln}P(i,j);$$

### Cys-Cys Interactions.

In principle, only 210 parameters are necessary for 20 amino acid residues. However, Cys-Cys contact requires special attention. Its propensity value is the largest compared to the other 209 contacts, because Cys residue tends to form disulfide bond with another Cys residue. Nevertheless, there are many Cys-Cys residue pairs that are in close spatial proximity but do not form disulfide bond. As a result, a mis-classification of a non-disulfide Cys-Cys contact as a disulfide bond will affect the overall score considerably, especially for small proteins with abundant Cys residues. The problem associated with the mis-classification of Cys-Cys contact was already pointed out in .

To avoid assigning the same score to the two different types of Cys-Cys contacts, we introduce a slightly more detailed propensity scores for Cys-Cys interactions. Since contacts between C:O, C:N, C:C, and O:O atoms never appear in disulfide bonds, a Cys-Cys contact pair lacking these atomic interactions is classified as a disulfide bond Cys-Cys contact if in addition the distance between their SG atoms is less than 2.5 Å. All other Cys-Cys interactions are classified as nonbonded Cys-Cys contact. The propensity values estimated for these two types of Cys-Cys contacts are listed in Table II.

### Bootstrap Resampling.

Because the sample size of 1,045 proteins in pdbselect  is limited, statistical modeling with approximations may be prone to errors. It is therefore essential to assess reliability of estimated contact potentials. Here we apply bootstrap techniques to calculate confidence intervals from simulated data sets . For alpha contact potential of a specific residue pair, (e.g., Trp-Trp), we denote the true value of the contact potential as $`\theta `$. Our probabilistic model (Equation 1) can be regarded as an estimator $`T`$ that gives the estimated value $`t`$ from the finite amount of data for $`\theta `$. Our goal is to calculate a $`95\%`$ confidence interval for $`\theta `$.

We resample 1,045 proteins independently $`R`$ times from the set of 1,045 proteins from Pdbselect , with replacement allowed. We have a simulated data set of $`Y_1^{},..,Y_R^{}`$, each contains 1,045 proteins. Some structures in the original Pdbselect  set appear multiple times, some appear once, and some never appear. We estimate the pair contact value $`\theta `$ from each of the $`R`$ samples, and obtain $`t_1^{},..,t_R^{}`$. For an equitailed $`95\%`$ confidence interval ($`95\%=12\alpha ,\alpha =2.5\%`$), we have the basic bootstrap confidence limits:

$$(t_{(R+1)(1\alpha )}^{},t_{(R+1)\alpha }^{})$$

The accuracy of these limits depend on $`R`$, and how well the distribution of $`\{t^{}t\}`$ agrees with that of $`T\theta `$. Perfect agreement occurs only when the distribution of $`T\theta `$ does not depend on any unknown variables.

To reduce possible errors due to unknown variables, we use studentized bootstrap. For the $`r`$th bootstrapped sample, we calculate:

$$z_r^{}=\frac{t_r^{}t}{v_r^{1/2}}$$

To obtain $`v_r^{}`$, we bootstrap with replacement again $`M`$ times the $`r`$th sample of the original bootstrap. We have:

$$v_r^{}=\frac{1}{M1}\underset{m=1}{\overset{M}{}}(t_m^{}\overline{t}^{})^2$$

where $`t_1^{},\mathrm{}t_M^{}`$ are calculated from the second bootstrap sampling. We then use the $`(R+1)\alpha `$th order statistic of the simulated values $`z_1^{},\mathrm{},z_R^{}`$, or $`z_{(R+1)\alpha }^{}`$ to obtain the studentized bootstrap confidence interval for $`\theta `$:

$$(tv^{1/2}z_{(R+1)(1\alpha )}^{},tv^{1/2}z_{(R+1)\alpha }^{})$$

Since $`M`$ bootstrap samples from the $`r`$th sample are needed to obtain $`v_r^{}`$, the total number of nested bootstrap samples is $`(M+1)R`$. We chose $`R=1,000`$ and $`M=50`$. Altogether, we generate $`1,001\times 50=50,050`$ bootstrap samples to calculate the confidence intervals for each of the $`209+2`$ pairwise alpha contact propensities.

### Database Selection.

In this study, we use Pdbselect  from http://www.cmbi.kun.nl /swift/pdbsel . Pdbselect  contains 1,045 proteins selected from the Protein Data Bank (PDB). The sequence identity between any pair of proteins in Pdbselect  is smaller than 25%.

## 3 Results

### Pairwise Alpha Contact Potentials.

The pairwise alpha contact propensities are listed in Table II. These propensities are calculated for all residue contacts that are at least 3 residues away in primary sequence. As expected, Cys-Cys has the highest propensities for contact interactions. Other residue pairs with the highest propensities for contact interactions ($`P(i,j)=1.42.5`$) are pairs of hydrophobic residues (e.g., Met-Met, Ala-Ala, Ile-Ile, Phe-Phe, Ile-Leu, and Ile-Met). The group of residue pairs with the second highest propensities ($`P(i,j)=1.21.3`$) are ionizable residues with opposite charges (e.g., Arg-Asp, Arg-Glu, Asp-Lys, and Gly-Lys). Residue pairs with lowest alpha contact propensities ($`P(i,j)=0.40.6`$) are dominated by pairs of ionizable residues of the same charge (e.g.,Arg-Lys, Asp-Glu, Lys-Lys, and Glu-Glu). The group of residue pairs with the second lowest alpha contact propensities ($`P(i,j)=0.50.7`$) are between ionizable residues and hydrophobic residues (e.g., Asp-Phe, Asp-Ile, Asp-Leu, and Glu-Val). Noticeably, pairs of Pro and ionizable/polar residues also have very low propensity for contacting interactions, probably due to the lack of a backbone-NH for H-bonding interactions.

The confidence intervals of these propensity values given by the studentized bootstrap procedure indicate that most of them are estimated accurately. Among the 209 parameters excluding Cys-Cys interaction, the 95% confidence intervals for 153 contact propensities are $`<0.2`$, a very tight interval. The confidence intervals of 36 contact propensities are $`<0.3`$. Contact propensities with the largest confidence intervals around 0.6 are Trp-Trp, Met-Met, Cys-Met, and Cys-Trp.

### Correlating and Clustering Similar Amino Acid Residues.

The overall behavior of pairwise contact interactions for a specific residue type is determined by its 20 pairwise contact propensity values. These values represent a profile of contact interactions specific for the residue type, and can be represented as a 20-dimensional vector $`𝒙`$.

We group the 20 types of amino acid residues by the Euclidean distance between the 20 vectors . Fig 2 shows the grouping of the 20 amino acid residues obtained by hierarchical clustering. As an exploratory tool for data analysis, hierarchical clustering can discover interesting and informative grouping patterns in data . In Fig 2, residues that have close contact propensity values to the 20 residue types are grouped together.

The pattern of groupings of residues clearly reflect the physical and geometric characteristics of the amino acid residues. Cys residue is different from all other residues, because of its propensity to form disulfide bond. The rest of the 19 residues can be broadly divided into two well defined branches of hydrophobic and hydrophilic residues. Among the hydrophilic residues, ionizable residues with positive charge and negative charge are grouped into two individual small clusters. Hydroxyl-containing residues (Ser and Thr) and amide residues (Asn and Gln) are also clustered into two individual small clusters. These residues are all capable of forming side chain hydrogen bonds, and their clusters are neighboring with each other, forming a larger cluster of Ser, Thr, Asn, and Gln. Among the hydrophobic residues, the branched residues (Val, Ile, and Leu), and small residues (Ala and Gly) are grouped into their own clusters. Aromatic residues Trp and Tyr also group into one cluster. Phe has strong hydrophobicity and is group with other strongly hydrophobic amino acid residues, rather than clustering with the other two aromatic residues. Pro and His are grouped together, probably because both do not form strong favorable contacts with either hydrophilic or hydrophobic residues.

The clustering pattern of residues by alpha contact propensity also resembles to a certain extent the clustering pattern derived from mutation matrix BLOSUM50  as reported by Murphy et al . For example, Arg with Lys, Asn with Gln, Ser with Thr, and Glu with Asp are all clustered tightly with each other. In addition to the distinct grouping of Cys in our clustering result, a notable difference is that by alpha contact propensity residues Ser, Thr and Pro are grouped with hydrophobic residues instead of hydrophilic residues.

### Discriminating Native Structures from Decoys.

An important method to determine the effectiveness of contact potential functions is to evaluate its success and failure in distinguish native protein structures from incorrectly folded decoy structures. We use three decoy data sets to assess alpha contact potential.

ProStar. The decoy sets in ProStar database  contain several subsets. The misfold subset contains conformations of 25 sequences which are obtained by placing these sequences on structures of different folds with the same number of residues. The conformations of the side chains are obtained by Monte Carlo sampling . Alpha contact potential succeeded in identifying all 25 native structures correctly (Table IIII).

For the Asilomar subset, alpha contact potential failed to identify 5 native structures out of 42 proteins. For the ifu subset, alpha contact potential failed for 20 out of 44 native structures. Alpha contact potential belongs to the class of residue-based potentials, similar to MJ, BT, CDF . As pointed out in , decoys in ifu subset are especially challenging for residue-based potentials, because they are conformations of short loops, and have only a small number of residue contact interactions.

In the subset pdberr and sgpa, the decoys are structures determined by diffraction but contain serious error, or are structures generated by molecular dynamics simulations starting from experimental conformations. Alpha contact potential missed one out of the five native structures.

Park & Levitt Set. This decoy test set contains native and near-native conformations of seven sequences, along with about 650 misfolded structures for each sequence. The positions of $`C_\alpha `$ of these decoys were generated by exhaustively enumerating ten selectively chosen residues in each protein using a 4-state off-lattice model. All other residues were assigned the phi/psi value based on the best fit of a 4-state model to the native chain. Conformations in the decoy sets all have low score by a variety of scoring functions, and have low RMSD to the native structure (Table IIIIII) .

The results of discrimination test are listed in Table IVIV and plotted in Fig 3. For five of the seven proteins, the native structures have lowest energy by alpha contact potential. For protein 3icb and 4rxn, the native structures have the 5th and 51st lowest energy values, respectively. For all proteins, decoys with the lowest energy are within 2.5 Å RMSD to the native structure.

The protein 3icb is a vitamin D-dependent calcium-binding protein. Although the energy of the native structure ranks as the 5th lowest energy, the RMSDs of the first four lowest energy decoys are all within 2.0 Å RMSD to the native structure. For a higher resolution structure (4icb at 1.60 Å) of this protein, the energy of 4icb by alpha contact potential is lower than any of the decoys. It is possible that this mis-classification might be due to the lower resolution of structure of 3icb.

Rubredoxin (4rxn) is an iron-sulfur protein. A Fe(III) ion is covalently bound in this structure with four Cys sulfur atoms, preventing these four Cys from forming two possible disulfide bonds. Because the protein description and the contact potential do not contain any information about the important covalent bonds of Cys with Fe-S cluster, it is reasonable to expect that native structure will not be of lowest energy if these important covalent bond interactions are unaccounted for. It is possible that the structure of rubredoxin might be different without the Fe-S cluster. The decoys at the lowest energy states form one or two fictitious disulfide bridges, and all are near-native structures with RMSD around 2 Å to the native structure. MJ and BT potential work better on 4rxn because they classify the four Cys-Cys contacts as disulfide bonds.

Lattice\_ssfit Set. The lattice\_ssfit set contains conformations for eight small proteins generated by ab initio protein structure prediction methods . The conformational space of a sequence was exhaustively enumerated on a tetrahedral lattice. A lattice-based scoring function was used to select the 10,000 best scoring conformations. Secondary structures were fitted to these conformations using a 4-state model . The 10,000 conformations were further scored using a combination of an all-atom scoring function , a hydrophobic compactness function, and a one-point per residue scoring function . The 2,000 best scoring conformations for each protein are selected as decoys for this data set.

The result (Table IVIV) shows that for six out of eight proteins, no decoy structures score better than the native structure. The exceptions are 1fca and 1trl. Similar to 4rxn in the Park & Levitt decoy set, ferrodoxin 1fca contains a Fe-S cluster. Its four Cys residues form four covalent bond with the four Fe(III) irons, instead of two disulfide bridges. These critical contacts again are unaccounted for in the alpha contact potential, and therefore the native structure of this protein was not identified successfully. 1trlA is a NMR solution structure of the C-terminal fragment (255-316) of thermolysin. NMR structures are far more difficult to recognize, as discussed in detail in . They are usually represented as an ensemble of conformations. The contact energies of conformations in the ensemble can be substantially different. It is conceivable that an energy function valid for crystal structures cannot reliably recognize native NMR structures . In addition, the structure 1trlA occurs in a dimeric state in the original PDB file. There is substantial interaction between the two chains. Because of this, it is unclear whether single subunit of 1trlA in monomeric state would retain the same conformation.

The decoy structures in this data set are generated by ab initio methods. None of them are near-native, and all have cRMSD to the native structure greater than 4.7 Å (Table IVIII). When decoys are so different from the native conformation, energy evaluated using alpha contact potential shows little correlation with the RMSD. The lowest energy decoys in this data set all have large RMSD, similar to results reported in  (data not shown).

Lmds Set. The local minima decoy set (lmds) contains decoys that were derived from the experimentally obtained secondary structures of ten small proteins that belong to diverse structural classes. Each decoy is a local minimum of a “hand made” energy function . Ten thousand initial conformations were generated for each protein by randomizing the torsion angles of the loop regions . The adjacent local minima were found by truncated Newton-Raphson minimization in torsion space. Each protein is represented in the decoy set by its 500 lowest energy local minima.

The alpha contact energy function works fairly well in the recognition of 1b0n-B, 1ctf, 1dtk, 2cro and 2ovo. 1dtk (Dendrotoxin K) contains three disulfide bonds in its native structure. However, in most of its 215 decoys, six Cys residues are spatially arranged together and form on average seven Cys-Cys contacts. For some decoys, up to 15 Cys-Cys contacts by distance cut-off can be found. The ability of discriminating disulfide bonded versus non-disulfide bonded Cys-Cys contacts probably makes the alpha contact potential discriminate better than the MJ and the BT potentials. The native structure of 1igd (IgG-binding domain from streptococcal protein G) ranks first by the MJ and the BT potential, but ranks 9th by alpha contact potential. The recognition of 1bba and 1fc2-C in this set failed for all residue-based contact potentials. 1bba is an atypical structure of a small protein determined by NMR, which forms a helix with random coil. 1fc2-C is a fragment of protein complexed with an immunoglobulin molecule. It possible that this protein may not maintain the same conformation without the complexed immunoglobulin.

By the criterion of the ranking of native protein, with the exception of the ion-sulfur proteins 4rxn and 1fca, the overall results shown in Table IVII and Table IVIV indicates that the performance of alpha contact potential in discriminating native protein from decoys is better than that of MJ and BT potentials for the misfold, ifu, asilomar, 4state\_reduced sets and the lattice\_ssfit set, and has comparable results for the lmds set, the pdberr set, and the spga set.

## 4 Discussion

### Contact Definition.

The alpha contacts introduced in this work are different from contacts by cut-off distances. Atoms in alpha contacts are all within a distance that depends on the identities of the two atoms. In this study, this distance is equal to the sum of the van der Waals radii of the two atoms, plus $`2\times 0.5`$ Å. Unlike contacts by distance cut-off, this distance is not a single fixed constant but depends on the atom types. Another important distinction of alpha contact is that only a subset of atoms satisfying the distance criterion will be counted as physical nearest neighbors, because we have an additional criterion: contacting atoms must have intersecting Voronoi cells. Alpha contacts represents the geometry more accurately and can capture contact interactions due to side chain size and orientation . In addition, no fictitious contacts will be introduced between two atoms when there is a third intervening atom . Perhaps this is the reason why alpha contact potential is sensitive to the presence of Fe-S clusters and other hetero atoms, which can be potentially exploited for determining whether a protein structure should contain hetero atoms.

For the 1,045 proteins in the pdbselect dataset, we compare contacts identified by distance cut-off with the threshold of two van der Waals atom radii plus $`2\times 0.5`$ Å and contacts identified by the alpha shape. We found that about 30 – 50 % of atom contacts detected by distance cutoffs are blocked by a third atom and hence do not have physical interactions. As a result, 3 – 6% residue contacts detected by distance cutoffs do not interact physically. Inclusion of these fictitious contact is problematic, especially in developing all-atom contact potentials, as well as in future studies when higher order interactions in the form of three or four body contacts are incorporated.

### Evaluating Discrimination of Alpha Contact Propensities by Bootstrap.

How robust is the results of decoy discrimination to the specific values of alpha contact potentials and the specific choices of the structures in the database? We further make use of the bootstrap resampling technique to evaluate the reliability of the discrimination results. As discussed earlier, we resampled 1,045 proteins in Pdbselect independently $`R=1,000`$ times with replacement allowed, and obtain 1,000 contact propensity matrices. Each is then used to discriminate the decoys in Table III. We use two parameters: $`\overline{r}`$, the average ranking of the native structure, and $`f`$, the times a native structure ranked as the structure with the lowest energy. Table IV shows that for many decoy sets (e.g., 1ctf, 1r69, 1sn3, 2cro, and 4pti in the 4state\_reduced set), the native structure always ranks first in the 1,000 bootstrapped energy evaluations. The performance using 1,000 different sets of “bootstrapped” potential values validates the robustness of the method deriving the alpha contact potential and the informativeness of the underlying protein structure database.

### Comparison with Contact by Geometric Centers.

Contact definition by distance cut-off is widely used in the development of many potential functions. Here we compare potentials obtained by alpha contact, denoted as $`H_\alpha `$, and by contact defined by cut-off distance between geometric centers of side chains, denoted as $`H_{gc}`$. Following reference , two residues are declared to be in contact if the distance $`d`$ between the geometric centers of their side chains is: 2Å$`<d<`$ 6.5Å. Geometric center based contact potential $`H_{gc}`$ are developed using the same Pdbselect data, with the same null model as that of alpha contact potential, and similarly counting only contact between residues that are three or more residues apart. Therefore, any difference between $`H_\alpha `$ and $`H_{gc}`$ is solely due to different geometric representation.

The $`\mathrm{log}`$ values of the parameters of $`H_\alpha `$ and $`H_{gc}`$ have an overall correlation coefficient of $`\rho =0.77`$. However, the contact maps of individual proteins by these two different contact definition are often substantially different. As an example, alpha shape dual simplicial complex gives significantly more contacts than cut-off distance by geometric centers of side chain for protein 1abe (Fig 4).

Fig 5 illustrates why such a discrepancy exists between these two contact models (see Figure 5 legend). The strong correlation between $`H_\alpha `$ and $`H_{gc}`$ is deceptive and this is reflected in another aspect. Although the correlation coefficient is high, the pairwise contact potentials may have very different values for $`H_\alpha `$ and $`H_{gc}`$. $`H_{gc}`$ categorically gives much higher propensity values for interactions between small residues. For example, $`H_{gc}`$ for Gly-Gly and Ala-Gly are 4.55 and 3.04, respectively, but $`H_\alpha `$ are only 1.48 and 1.39, respectively. Gly-Pro interaction is strongly favorable by $`H_{gc}`$ ($`P(Gly,Pro)=1.84`$), but is unfavorable by $`H_\alpha `$ (0.87). On the other hand, $`H_{gc}`$ gives much lower propensity values for interactions between large residues. For example, $`H_{gc}`$ for Trp-Trp and Phe-Trp are 0.39 and 0.56, respectively, but $`H_\alpha `$ are 1.75 and 1.65, respectively. In addition, many pair contact interactions between Trp and another residue with large side chains (such as Arg, Tyr, Phe, His, Leu, Ile, and Met) are unfavorable ($`<1.0`$) by $`H_{gc}`$, but are favorable by $`H_\alpha `$.

These differences lead to different discrimination in identifying native and near-native protein structures from decoy structures. Because it is impossible to define refined potential for Cys-Cys contacts in $`H_{gc}`$ model as in alpha contact model, we exclude proteins containing Cys-Cys contacts to avoid complication for comparison. Table IVIV shows that $`H_{gc}`$ can only recognize two native structures out of nine proteins in the union set of the 4-state decoy set and the lattice-ssfit decoy set. In addition, the $`z`$ scores for native structures are generally higher when using $`H_\alpha `$ than $`H_{gc}`$. For the 4-state decoy set, the correlation coefficient $`\rho `$ by $`H_\alpha `$ is also higher. These results indicate that geometric description of protein structures are important, and contact model by alpha shape are more accurate with more discriminative information for identifying native-like structures. However, the cut-off distance approach may be more convenient to implement, and it is possible to gain further improvement by setting the values of the cut-off threshold as a variable depending on the type of contacting residue pairs.

### Comparison with Other Potentials.

We further compare alpha contact potential with several previously developed residue contact potentials, including MJ , BT , SK , and TD  potentials. We take the values of MJ potential from Table 3 of reference . Following the author’s recommendation, the average hydrophobicity ($`e_{rr}`$ = -2.55) was subtracted from the potential, and a pair of residues is declared to be in contact if the geometric centers of their side chain is within an interval of 2.0 Å to 6.5 Å. The values of BT potential is taken from  and was obtained by rescaling MJ potentials with Thr as the reference solvent . The correlation coefficients $`\rho `$ and dispersions  between each pair of potentials are shown in table Table VV. The correlation coefficients $`\rho `$ for alpha contact potential $`\mathrm{log}P(i,j)`$ and these residue contact potentials are: $`\rho =0.66,0.80,0.61`$, and $`0.66`$ for MJ, BT, SK, and TD potentials, respectively. The dispersion as defined in  (p.363, Formula 4) are $`1.45,0.28,0.51`$, and $`0.39`$, respectively. Because the contact map obtained by alpha edges can be substantially different from other contact definition (Fig 4), the absolute value of energy by alpha contact potential and by other potentials for the same structure can be substantially different.

### Long Range Interactions.

Interactions between residues with large sequence distance $`d`$ are relatively rare. We found that they occur more likely in the interior of a protein than on the surface (data not shown). Identifying such interactions are of particular interest, because these interactions result in significant reduction of conformational entropy. Prediction of protein structures seem to be most difficult for those with large contact order , namely, those with significant amount of interactions between residues with large sequence distances.

The bootstrap procedure introduced here provides a method to reliably identify the contact pairs of long sequence separation whose propensity values can be confidently assessed (see Supplementary Online Material for tables of alpha contact potential for $`d30`$). Among all possible 210 pairwise interactions, 9 contact pairs with high propensity (lower value of 95% confidence interval $`>1.5`$) can be reliably assessed for $`d>30`$. In addition to Cys-Cys, they include hydrophobic-hydrophobic interactions (Gly-Gly, Met-Met, Ile-Ile, Phe-Phe, Val-Val, Met-Phe), salt-bridge interactions (Arg-Asp, Asp-Lys), and Pro-Trp.

Some long range interactions are clearly associated with specific secondary structures. After correction for prior probability to be in a particular secondary structure, we found that Met-Met contact has a high propensity to occur between two helices (h) or two beta-strands (s) ($`P(Met,Met)_{hh}=2.4,P(Met,Met)_{ss}=2.3`$), and a low propensity to occur between either a helix and a coil (c) ($`P(Met,Met)_{hc}=0.65`$), or a strand and a coil ($`P(Met,Met)_{sc}=0.56`$). Similarly, because Gly is a helix breaker, long-range Gly-Gly contact has a high propensity to occur between two coils ($`P(Gly,Gly)_{cc}=3.2`$), and low propensity to occur between two helices ($`P(Gly,Gly)_{hh}=0.53`$).

### Reduced Alphabet for Amino Acid Residues.

The clustering of pairwise alpha contact potentials shown in Figure 2 suggests that many residues behave similarly in contact interactions. This points to possible degeneracy of the amino acid alphabet . Reduced alphabet is important because a smaller size of alphabet will lead to exponentially more efficient sampling methods in sequence design and protein engineering . Many experimental and computational studies already suggested that a minimum number of amino acid residue types far less than 20 may be adequate for protein folding . Wang and Wang examined different ways to reduced MJ interaction matrix and concluded that by minimizing mismatches, a reduced alphabet of just five amino acid residue types can be used to construct sequences with good foldability and kinetic accessibility . The reduced 5 alphabet set coincide with the same alphabet set reported in the work of Riddle et al, where fully functional constructs for a small 57-residue beta-barrel protein can be experimentally obtained when residues in 38 out of 40 selected amino residues are drawn from the alphabet set of I, K, E, A, and G. Murphy et al further examined reduced alphabets based on BLOSUM50 substitution matrix . When using a variety of reduced alphabets with size ranging from 10 to 20, they found that there is little loss of the information necessary to select structural homologs in a database of representative protein sequences using dynamic programming based global alignment.

We continue investigation in this direction and study the capability of various reduced alphabet sets in discriminating native proteins from decoys. Fig 2 provides a natural way of reducing the residue alphabet set that is similar to the approach used by Murphy et al . By placing a horizontal line at different vertical heights, we can obtain a reduced residue alphabet that is determined by the heights of the branching points in the dendrogram from the hierarchical clustering. For example, we can have an alphabet of 7 residue types at height about 1.5: $`𝒜`$ = {D,E}, $``$ = {R,K}, $`𝒞`$= {S,T, N, Q, H, P}, $`𝒟`$ = {V, I, L, M, F, W}, $``$ = {W, Y}, $``$ = {A, G}, and $`𝒢`$ = {C}. An alphabet of two residue types would take Cys as a residue type, and everything else grouped into another residue type. An alphabet of three residue types would have Cys, polar residues, and hydrophobic residues.

Does reduced alphabet still capture the basic information of protein contact interactions? We use Equation 1 to estimate pairwise alpha contact potentials for reduced alphabet size of 2, 3, 5, 7, 8, 9, 10, 11, 15, and 20, and test their effectiveness in selecting native structure from decoys in the Park and Levitt data set. Fig 6 shows the result when an alphabet set of 9 residue types, plus the three types of Cys-Cys contacts are used. Remarkably, the discrimination of native conformation from decoys is almost as good as when there are 20 different residue types.

Figure 7 shows the detailed results of discriminating native structure of 3icb from decoys using potentials derived from different alphabet sets. The average RMSDs of $`n`$ decoys with lowest energy by potentials of different alphabets are calculated. Smaller average RMSD of these lowest energy decoys indicates that a large fraction of them are near-native structures. This would suggest good discrimination. The top $`n=100`$ decoys of lowest energy are all found to have very similar average RMSD to the native structure, regardless the size of the alphabet. This suggests that alphabet with just a few residues have respectable results in decoy discrimination. Table VIVI further shows the correlation coefficient of energy and RMSD for all proteins in the Park and Levitt data set when using different alphabets. The results indicate that an alphabet of 7 residues would have very similar performance as an alphabet with 20 residues in discriminating decoys. Our results extended earlier work where subjectively defined alphabet sets were used to extract contact residue potentials by an iterative optimization method . It was found that potential derived using such reduced alphabets were effective in discriminating decoys generated by gapless threading. Here we showed that similar conclusion can be drawn for statistical potential using alphabet sets derived from natural clustering of residues, in discriminating natives against more stringent compact decoy conformations generated by off-lattice model . Our conclusion is also consistent with a recent study where it is shown that much of the information in pairwise contact potential is related to just a few variables such as hydropathy, charge, disulfide bonding, and residue burial .

Although the alphabets we used have different number of residues, they are all developed with one aspect in common, i.e., the contacts are all derived from the dual simplicial complexes, which provides a faithful representation of the geometry. This suggests that as long as the same space filling pattern is conserved, the specific residue types are not critical in many cases. It seems that packing geometry plays a very important role, but the specific residue types are often replaceable. This observation is consistent with experimental results where it is well known that proteins are robust against many mutations.

### Edge Simplices and Tetrahedron Simplices.

Pairwise alpha contact potential only considers the edge simplices or 1-simplices in the dual simplicial complex. There have been several studies of statistical potential based on 3-simplices or tetrahedra . In the work of Tropsha et al, 3-simplices are obtained from unweighted Delaunay triangulation . In these studies, all residues are treated as balls of equal size located at Cα or Cβ positions, and a cut-off distance is used to remove tetrahedra that are considered to be too large. Our approach is different. Our model is instead based on the weighted dual simplicial complex, a different simplicial complex formed by a subset of the simplices from the weighted Delaunay triangulation of all atoms in the molecule. The dual simplicial complex or the alpha shape allows modeling at atomic level. Therefore, in our approach contacts can be defined by the full side chain and main chain atoms. Additionally, atoms are assigned with appropriate non-uniform van der Waals radii . Finally, because the dual simplicial complex reflects the precise contact geometry, we avoid the use of heuristic cut-off thresholds necessary to eliminate a subset of the simplices from the Delaunay triangulation. Our contacts represents accurately geometry of the structure. We discussed earlier the differences between the alpha contact and contact by distance cut-off between geometric centers of side chains. It is conceivable that similar difference will result between alpha contact and the approach described in .

### Summary.

In this work, we introduced a novel representation of protein structures using edge simplices of the alpha shape, or the dual simplicial complex of the protein structure. By describing pairwise contact interactions with simplicial edges, we developed alpha contact potential based on the statistics of edge simplices. We also developed a bootstrap model which provides confidence interval estimations, including those of long range interactions. We found that alpha contact potential performs well in decoy structure discrimination. By comparing with alternative contact potential, we conclude that geometric representation of contact interaction is important, but the specific residue types are often interchangeable.

## 5 Acknowledgment

We thank Dr. Hui Lu for very helpful discussion. This work is supported by funding from National Science Foundation CAREER DBI0133856, DBI0078270 and MCB998008 and American Chemical Society/Petroleum Research Fund 35616-G7.


format: text

text:
# Electronic structure of half-metallic double perovskites

## I Introduction

The interest in double perovskites A2B’B”O6, with A being alkaline earth (AE) ion, and the transition metal sites (B-sites) alternately occupied by different cations B’ and B”, has really taken off with the paper by Kobayashi et al. kobayashi , demonstrating a convincing half-metallic behaviour in Sr2FeMoO6. The calculated electronic structure of this compound shows a gap in the majority spin component, while for the minority spin component one observes strongly hybridized Fe 3$`d(t_{2g})`$, Mo 4$`d(t_{2g})`$, and O $`2p`$ states at the Fermi level. Their study has also established that the latter compound exhibits tunnelling magnetoresistance (TMR) at low applied magnetic field even at room temperature. The double perovskites were discovered in 1960s longo ; sleight ; patter ; galasso ; sleight1  and seem to be metallic and ferromagnetic for B’=Fe, and B”=Mo or Re. Unlike the manganites, they are stoichiometric and have large Curie temperatures, Tc, which makes them very promising for technological applications. The Sr2FeMoO6 compound has Tc of about 450 K, but also Ba2FeMoO6 and Ca2FeMoO6 have been shown to have interesting TMR characteristics and high Tc’s.coey  Although Sr2FeReO6 (see Ref. SrRe, ) shows similar TMR properties to those of Sr2FeMoO6, recent studies claim that the Ca-compound, Ca2FeReO6, unlike Ca2FeMoO6, may be insulating.Millis  It has been suggested that the magnetic temperature Tc, and whether the ground state is metallic or insulating, vary as A is changed from Ba to Sr to Ca. Millis1  Also mis-site disorder may have pronounced effect on the magnetic properties of these compounds. ogale ; Tanusri

Recent experimentsbesse  have established in Sr2FeMoO6 an electronic structure and configuration of five localized majority electrons (high spin state) on Fe and one delocalized electron shared between Mo and the other sites. We aim to validate this electronic structure by a more systematic study involving two other alkaline earths, Ca and Ba, as well as Sr. The aspect of a delocalized Mo electron is investigated by replacing Mo by Re and investigating the amount of delocalized electrons shared between Re and the other sites.

In this paper we present an application of the self-interaction corrected local spin density (SIC-LSD) approximation to the double perovskites Ba2FeMoO6, Ca2FeMoO6, Sr2FeMoO6, and Ca2FeReO6. Specifically, we concentrate on the electronic and magnetic properties of these compounds and in particular the size and relative orientation of the spin moments of Fe and Mo (or Re). In all the compounds the SIC-LSD calculations find a spin moment of about -0.4 $`\mu _B`$ on the Mo and about -1.1 $`\mu _B`$ on Re, with an opposite orientation to the respective spin moments on Fe sites. In addition, we find these double perovskites to be half-metallic, with a well defined gap in the spin-up density of states, and strong hybridization at the Fermi energy between the spin-down Fe 3$`d`$, Mo 4$`d`$ (or Re 5$`d`$), and O 2$`p`$ states.

The earlier first-principles calculations performed for double perovskites had to be done with GGA and implement theoretical optimisation of oxygen positions to obtain the half-metallic groundstate. kobayashi ; sarma  Other applications were done with the LDA+U, fang  which similarly to SIC-LSD provides better description of localized states in solids. The present study is the first SIC-LSD application to double perovskites. Its advantage is that it treats both localized and delocalized states on equal footing, and therefore should be well suited for describing the localized nature of Fe $`3d`$ electrons as well as the Fe valency.

The remainder of the paper is organized as follows. In the next section, we briefly present the theoretical background of the SIC-LSD approach. In section III, the technical and computational details regarding the application to double perovskites are elaborated upon. In section IV, we present total energy and density of states calculations. There we also compare magnetic moments and determine valency of Fe in the studied compounds. The paper is concluded in section V.

## II Theory

The basis of the SIC-LSD formalism is a self-interaction free total energy functional, $`E^{SIC}`$, obtained by subtracting from the LSD total energy functional, $`E^{LSD}`$, a spurious self-interaction of each occupied electron state $`\psi _\alpha `$PZ81 , namely

$$E^{SIC}=E^{LSD}\underset{\alpha }{\overset{occ.}{}}\delta _\alpha ^{SIC}.$$
(1)
Here $`\alpha `$ numbers the occupied states and the self-interaction correction for the state $`\alpha `$ is

$$\delta _\alpha ^{SIC}=U[n_\alpha ]+E_{xc}^{LSD}[\overline{n}_\alpha ],$$
(2)
with $`U[n_\alpha ]`$ being the Hartree energy and $`E_{xc}^{LSD}[\overline{n}_\alpha ]`$ the LSD exchange-correlation energy for the corresponding charge density $`n_\alpha `$ and spin density $`\overline{n}_\alpha `$. It is the LSD approximation to the exact exchange-correlation energy functional which gives rise to the spurious self-interaction. The exact exchange-correlation energy $`E_{xc}`$ has the property that for any single electron spin density, $`\overline{n}_\alpha `$, it cancels exactly the Hartree energy, namely

$`U[n_\alpha ]+E_{xc}[\overline{n}_\alpha ]=0.`$ (3)
In the LSD approximation this cancellation does not take place, and for truly localized states the correction (Eq. (2)) can be substantially different from zero. For extended states in periodic solids the self-interaction vanishes. Consequently, the SIC-LSD approach can be viewed as a genuine extension of LSD in the sense that the self-interaction correction is only finite for spatially localized states, while for Bloch-like single-particle states $`E^{SIC}`$ is equal to $`E^{LSD}`$. Thus, the LSD minimum is also a local minimum of $`E^{SIC}`$.

It follows from minimization of Eq. (1) that within the SIC-LSD approach such localized electrons move in a different potential than the delocalized valence electrons which respond to the effective LSD potential. For example, in the case of the double perovskites, five (Fe3+) or six (Fe2+) Fe $`d`$ electrons move in the SIC potential, while all other electrons feel only the effective LSD potential. Thus, by including an explicit energy contribution for an electron to localize, the ab-initio SIC-LSD describes both localized and delocalized electrons on an equal footing, leading to a greatly improved description of static Coulomb correlation effects over the LSD approximation. Assuming various atomic configurations, consisting of different numbers of localized states, one can explore the corresponding local minima of the SIC-LSD energy functional of Eq. (1), and determine the lowest energy solution and valency. The advantage of the SIC-LSD formalism is that for such systems as transition metal oxides or rare earth compounds the lowest energy solution will describe the situation where some single-electron states may not be of Bloch-like form. Regarding the valency, it is determined here as the integer number of electrons available for band formation, i.e.,

$$N_{val}=ZN_{core}N_{SIC},$$

where $`Z`$ is the atomic number (26 for Fe), $`N_{core}`$ is the number of core (and semi-core) electrons (18 for Fe), and $`N_{SIC}`$ is the number of localized, i.e., self-interaction corrected, states (either five or six, respectively for Fe3+ and Fe2+).

In the present work the SIC-LSD approach has been implementedTSSW98  within the linear muffin-tin-orbital (LMTO) band structure methodoka75  in the tight-binding representationAJ84 , where the electron wave functions are expanded in terms of the screened muffin-tin orbitals, and the minimization of $`E^{SIC}`$ becomes a non-linear problem in the expansion coefficients. The atomic spheres approximation (ASA) has been used, according to which the polyhedral Wigner Seitz cell is approximated by slightly overlapping atom centered spheres, with a total volume equal to the actual crystal volume.

## III Calculational details

All the compounds studied here crystallize in ordered double perovskite structures. The simplest system is Ba2FeMoO6 which occurs in the cubic structure, with the Fm3m space group. The Sr2FeMoO6 compound has the body-centered tetragonal bct structure, with the I4/mmm space group. Replacing Sr by a smaller still Ca element results in a monoclinic structure for the Ca2FeMoO6 and Ca2FeReO6 compounds, with the P21/n space group. The experimental lattice parameters used in the calculations for the first three compounds have been taken from Refs. coey,  and kobayashi, , and for the last compound from Ref. Millis, . For the Ca compounds two formula units (20 atoms) have been necessary. In addition, for a better space filling also eight empty spheres have been introduced in Re-compound. Concerning the linear muffin-tin basis functions, we have used $`4s`$, $`4p`$, and $`3d`$ partial waves on all Fe atoms, and treated them as low-waves, on the oxygens ($`2s`$ and $`2p`$) and Ca ($`4s`$ and $`3p`$) in addition to the $`s`$ and $`p`$ low partial waves, $`3d`$partial waves have been treated as intermediate. On the empty spheres the $`1s`$waves have been considered as low- and the $`2p`$’s as intermediate-waves. On the Ba ($`6s`$, $`5p`$, and $`5d`$), Sr ($`5s`$, $`4p`$, and $`4d`$), Mo ($`5s`$, $`5p`$, and $`4d`$) and Re ($`6s`$, $`6p`$, and $`5d`$) atoms, the $`s`$, $`p`$, and $`d`$waves have been treated as low and the $`4f`$waves as intermediate. In all the cases, the SIC-LSD calculations have been started using the self-consistent LSD charge densities with the Perdew and Zunger exchange-correlation potential. PZ81  Mostly the Fe 3$`d`$ electrons have been treated as localized states (in some test calculations we have also treated Mo $`d`$’s as localized), meaning that they have been moving in the self-interaction corrected LSD potentials, while all other electrons have seen the effective LSD potential. In all the cases studied, we have allowed for a number of different configurations of SIC 3$`d`$ states to find the ground state solution. Namely, treating five 3$`d`$ electrons on the Fe’s as localized, moving in the SIC-LSD effective potentials, meant that the Fe ions were trivalent. Similarly, for divalent Fe’s, six 3$`d`$ electrons on the Fe’s were considered to be localized. By comparing total energies corresponding to different localized Fe 3$`d`$ configurations we have been able to find the ground state configuration which for all the compounds has been that of Fe3+. However, other configurations have been energetically close, with the Fe2+ being insulating in all Mo-based compounds. Although, unlike in the work by Kobayashi et al., kobayashi  optimizing oxygen positions in Sr2FeMoO6 was not needed for obtaining half-metallic ground state, however in order to study its effect on the magnetic properties and to make contact with Ref. kobayashi, , for this compound we have also performed calculations for the scenario where the theoretically optimized positions of Ref. kobayashi,  have been implemented. besse

## IV Results and discussion

The SIC-LSD spin-polarized total densities of states (DOS) per formula unit are presented in Fig. 1 for Ca2FeReO6, Ca2FeMoO6, Sr2FeMoO6, and Ba2FeMoO6 compounds. There we also show the O $`2p`$ spin polarized densities of states. The spin polarized densities of states of Fe $`3d`$ and Mo $`4d`$ and Re $`5d`$ are given in Fig. 2, for all the different compounds. First thing to note is the similarity between the densities of states, independently of which alkaline earth compound is considered. For all compounds we see convincing half-metallic behaviour, with well defined gap at the Fermi energy in the majority spin channel, and strongly hybridized Fe $`3d`$, Mo $`4d`$ (or Re $`5d`$) and oxygen $`2p`$ states, in the other spin channel. The half-metallic character is reflected in the total spin magnetic moments that are all integer, as can be seen in Table I, where we summarize the total and species decomposed spin moments for all the compounds together with the respective volumes per formula unit. As can be seen in Fig. 1, the valence band consists predominantly of O $`2p`$ states, but contains a small admixture of Mo, Re- and Fe $`d`$states. The unoccupied majority Re $`d`$ bands lie closer to the Fermi energy in comparison with the unoccupied majority Mo $`d`$ states in Ca2FeMoO6. We can also see that the unoccupied Mo $`d`$ states move closer to the Fermi energy as one traverses the series in the sequence Ca, Sr and Ba. However, the majority spin channel maintains a band gap for all these compounds. In the minority spin channel these systems are metallic with O 2p, Mo 4d and Fe 3d states straddling the Fermi energy. These are the states which reduce the Fe magnetic moment from essentially 5 $`\mu _B`$ to 3.87 $`\mu _B`$ (Ca2FeReO6), 3.76 $`\mu _B`$ (Ca2FeMoO6), 3.71 $`\mu _B`$ (Sr2FeMoO6) and 3.81 $`\mu _B`$ (Ba2FeMoO6) (see Table I). From these DOS one can also expect the oxygen to acquire a moment. This happens according to Table I and it is parallel to the Fe spin moment and is non-neglegible at around 0.1 $`\mu _B`$. However it is the Mo and Re sites which are of interest. Substantial spin magnetic moments of -0.4 $`\mu _B`$ and -1.1 $`\mu _B`$ are found on the Mo and Re sites respectively. This results from the strong hybridization with Fe in the minority spin channel. This also explains the antiparallel alignment of these spin magnetic moments with respect to Fe.

The densities of states in the majority spin channels for these four compounds are insulating, since a fully occupied O $`p`$ band and five localized Fe $`d`$ states accomodate 23 (6$`\times `$3+5=23) valence electrons. Of the remaining 19 or 20 valence electrons for the Mo and Re compounds respectively, 18 fill up the minority O $`p`$ band. This leaves one or two electrons for the Mo and Re compounds to occupy the minority states which straddle the Fermi level, and have Fe $`d`$, Mo or Re $`d`$ and oxygen $`p`$ character. We have also considered a scenario of localising through the self-interaction correction one Fe (Mo) $`d`$ state in the minority spin channel describing Fe$`{}_{}{}^{2}+`$ (Mo$`{}_{}{}^{5}+`$) respectively. In this case the states straddling the Fermi energy which accomodate this one electron are pulled down below the bottom of the O $`p`$ valence bands and we also obtain an insulating state in the minority spin channel. However, these scenarios are energetically unfavourable by 70 mRy and 103 mRy for Mo5+ and Fe2+ respectively. Obviously the gain in localisation energy does not compensate the loss in hybridization energy.

It is interesting to see how little the magnetic properties of these compounds are affected by the size of the alkaline earth atom. There is hardly any variation in the size of the Fe and Mo or Re moments with the substantial change in volume of the compounds (see Table I). Also the reduction in the size of the Fe and Mo spin moments for Sr2FeMoO6, as a result of implementing theoretical optimization of oxygen positions, as given by Kobayashi et al.,kobayashi  does not seem too important although brings the size of our calculated moments close to the values of Ref. kobayashi, . The reduced values are given in parenthesis in the fourth column of Table I. Concerning spin moments of the other species, the small changes cannot be resolved when rounded off to two digits after the decimal point.

Finally, we would like to comment on the valency of Fe and the ground state configuration for all the compounds. Regarding energetics, in all SIC-LSD calculations the Fe3+ configuration has been most energetically favourable, followed by Fe4+ (energetically unfavourable by 53, 61 and 83 mRy in Ba-, Sr- and Ca-compounds, respectively), and Fe2+ solution, at least twice as unfavourable as the former configuration (specifically, by 103, 122 and 250 mRy). Moreover, treating also Mo $`d`$ states as localized has been energetically unfavourable, and the more states have been localized the more unfavourable solutions have been obtained.

## V Conclusions

In summary, the SIC-LSD approach has been applied to four different double perovskites, and in all cases the half-metallic ground state solutions have been obtained. This seems to be a rather generic result for double perovskites involving Fe with a filled majority d shell and transition metal ions which donate few, one (Mo) to two (Re), d electrons. This type of materials could therefore be a fruitful playground for the discovery of new spintronics materials. In particular, also the Ca2FeReO6 has been found to be half-metallic, and not insulating as indicated in Ref. Millis, . The electronic structure for all studied compounds appears to be very similar, with a gap in the majority spin component, and strongly hybridized bands at the Fermi level for the minority spin component. The spin magnetic moments induced by hybridization on the Mo sites have been of the similar magnitude, about -0.4 $`\mu _B`$, independently of the alkaline earth element, while that of Re is of the order of -1.1 $`\mu _B`$. These induced spin magnetic moments are antiparallel aligned with the Fe spin moment. Concerning the Fe valency, for all the compounds the trivalent configuration has been most favourable, followed by the tetravalent and divalent ones.

## Acknowledgements

We thank Dr. Tanusri Saha-Dasgupta for useful communications and Professor P.H. Dederichs for valuable discussions.


text:
# On Stellar Dynamo Processes and Differential Rotation

## 1 Introduction

During the last three decades considerable evidence for stellar magnetic activity has been gathered (Wilson 1978, Noyes et al. 1984, Baliunas et al. 1995). These observations based on CaII emissions, revealed that some stars exhibit irregular magnetic activity while others have solar-like activity cycles. What could be the source of such magnetic fields? Are the different types of activity (cyclic vs. irregular) indications of differences in the generation and amplification of magnetic field? Based on the ohmic diffusion time $`\tau _{ohm}R^2/\eta `$ (where $`R`$ is the stellar radius and $`\eta `$ is the magnetic diffusivity, leading to $`\tau _{ohm}10^{10}`$ years for the Sun), many of those stars could have retained a fossil magnetic field over their main sequence lifetime without running a dynamo (Weiss 1994). However, it is tempting to search for a link between stellar global properties such as mass, rotation or convection and stellar magnetic activity. It is indeed important to notice that convection is always present in stars and that most of them rotate. Stars less massive than $`1.3`$ solar mass $`M_{}`$ display convective envelope whereas the more massive ones possess a convective core. The latter appear due to an increasing contribution to the star’s energy generation of the CNO nuclear cycles. Further, the range of rotation rates varies significantly from very slow rotators ($`P>150`$ days) to near breakup velocities. One general trend is that older stars rotate slower. Another is that the presence of magnetic fields increases the probability for a star to be a slow rotator. More precisely Durney & Latour (1978) have shown that stellar magnetic activity for main sequence stars can reasonably been explained by evaluating the Rossby number $`R_o`$. In their model, strong magnetic activity is seen in stars for which $`R_o<<1`$ (i.e. their convective motions strongly feel rotation). In reality stars are more complex, and there is evidence that slowly rotating Ap stars exhibit stronger magnetic field than fast rotating A stars (Moss 2002). Here even though A stars have a convective core that is likely to generate a dynamo field, some quenching mechanism has to be invoked to reproduce the observations–unless the observed field in Ap star is of fossil origin and the dynamo generated field screened by the radiative envelope. I will here focus my attention on solar type stars (see Mestel 2001 for stars of other spectral type). In particular, I will examine the dynamical establishment of the solar differential rotation and meridional circulations and their nonlinear link with magnetic fields. I will present recent results from simulations of solar convection within full 3–D spherical shells, discussing the angular velocity $`\mathrm{\Omega }`$ profiles that can be achieved in the bulk of such convection zones and the level of dynamo induced magnetism that can be sustained there.

## 2 Our Numerical Approach

The anelastic spherical harmonics (ASH) code solves the 3–D MHD anelastic equations of motion in a rotating spherical shell geometry using a pseudo-spectral/semi-implicit method (Miesch et al. 2000, Brun & Toomre 2002,2003). The model is a simplified description of the solar convection zone: solar values are taken for the heat flux, rotation rate, mass and radius, and a perfect gas is assumed since the upper boundary of the shell lies below the H and He ionization zone. The computational domain extends from 0.72 $`R_{}`$ to 0.97 $`R_{}`$, thereby concentrating on the bulk of the unstable zone and here not dealing with penetration into the radiative interior. The effects of the steep entropy gradient close to the surface has been softened by introducing a subgrid scale (SGS) transport of heat to account for the unresolved motions, and enhanced eddy diffusivities are used in these large eddy simulations (LES). The typical density difference across the shell in radius is about 30.

## 3 Solar Convection and Rotation

### 3.1 Convective Patterns

Figure 1 displays the evolution of enstrophy (vorticity2) in the case called $`E`$ (Brun & Toomre 2003) over 10 days in time near the top of the domain. The vantage point is in the uniformly rotating frame used in our simulations. The convection patterns are highly time dependent. Some of the pattern evolution is related to the advection by zonal flows associated with the differential rotation driven by the convection relative to this frame (i.e fast equator/slow poles). There is an asymmetry between the broad upflows showing a weaker vortical content and the narrow fast downflows at their periphery. This leads to a downward transport of kinetic energy. The strong correlations between warm upward motions and cool downward motions are essential in transporting the heat outward.

Vortical structures are found at the interstices of the downflow network, although strong fronts are seen in a band near the equator. The strongest of these vortex tubes or ‘plumes’ extend through the whole domain depth. These plumes represent coherent structures that are surrounded by more chaotic flows. They are counterclockwise in the northern hemisphere and clockwise in the southern one. They tend to align with the rotation axis and to be tilted away from the meridional planes, leading to Reynolds stresses that are crucial ingredients in redistributing the angular momentum within the shell (Brun & Toomre 2002).

### 3.2 Differential Rotation

The differential rotation profile in latitude and radius associated with the vigorous convection of case $`E`$ is shown in Figure 2. For simplicity, we have converted the mean longitudinal velocity $`\widehat{v}_\varphi `$ into a sidereal angular velocity $`\mathrm{\Omega }`$, using $`\mathrm{\Omega }_o/2\pi =414`$ nHz (or 28 days) as the reference frame rotation rate. In the contour plot, the near polar regions have been omitted due to the difficulty of forming stable averages there, since the averaging domain is small but the temporal variations large. Case $`E`$ exhibits a fast (prograde) equatorial region and slow (retrograde) high latitude regions. This is due to correlations in the velocity components leading to significant Reynolds stresses. These Reynolds stresses are intimately linked to the influence of Coriolis forces acting upon the convecting motions and to the presence of plumes tilted both away from the local radial direction and out of the meridional plane. Such correlations have been identified in local high resolution Cartesian domains as well (Brummell et al. 1998). These lead to an equatorward transport of angular momentum, resulting in the slowing down of the high latitude regions and speeding up of the equatorial zone. At low latitudes there is some alignment of $`\mathrm{\Omega }`$ contours with the rotation axis. At mid latitudes, the angular velocity is nearly constant along radial lines, in good agreement with helioseismic deductions (Schou et al 1998). Further, case $`E`$ exhibits a monotonic decrease of $`\mathrm{\Omega }`$ with increasing latitude, a property that has been difficult to achieve in 3–D spherical convection calculations. Indeed, most other cases have their equator to pole contrast $`\mathrm{\Delta }\mathrm{\Omega }`$ confined to mid latitudes extending to about $`42^{}`$ (Brun & Toomre 2002). The differential rotation contrast between the equator and $`60^{}`$ in case $`E`$ is 110 nHz (or 26% relative to the frame of reference), thus being very close to the 92 nHz (or 22%) variation observed in the Sun. A sense of the asymmetry present in case $`E`$ can be assessed both in the contour plot and in the latitudinal cuts (right panel of Fig. 2), where we have plotted $`\mathrm{\Omega }`$ in the north (dotted) and south (dashed) hemispheres, along with their mean. The convection itself exhibits some asymmetry between the two hemispheres (cf. Fig. 1), and so it is not surprising that the mean flows driven by the convection do the same. These asymmetries are expected to diminish over longer temporal averages. Mean field models of the solar differential rotation (Kichatinov & Rüdiger 1995) have advocated that a thermal wind balance (involving pole to equator temperature contrasts) could be the cause of the non-cylindrical profile in $`\mathrm{\Omega }`$. This could come about through the baroclinic nature of the convecting motions yielding some latitudinal heat fluxes, resulting in the breakdown of the Taylor-Proudman theorem. Although it is indeed true that case $`E`$ exhibits latitudinal variation of entropy and temperature fluctuations relative to the mean, these are not the most dominant players everywhere in the shell. A temperature contrast of few degree K seems compatible with a $`\mathrm{\Delta }\mathrm{\Omega }/\mathrm{\Omega }_o`$ of $`30\%`$. However, we find that the Reynolds stresses are the main agents responsible for the equatorial acceleration achieved in our simulations, and thus the solar differential rotation is dynamical in origin.

### 3.3 Meridional Circulation

The meridional circulation associated with the vigorous convection in case $`E`$ is maintained variously by Coriolis forces acting on the differential rotation, by buoyancy forces, by Reynolds stresses and by pressure gradients, and thus can be thought as departures from a simple geostrophic balance (Brun & Toomre 2002). The meridional circulation seen in case $`E`$ exhibits multi-cell structure both in latitude and radius, and given the competing processes for its origin, it is not straightforward to predict. Typical amplitudes for the velocity are of order 25 m/s, comparable to local helioseismic deductions (Haber et al. 2002). The flow is directed poleward at low latitudes, with return flow deeper down. The temporal fluctuations in the meridional circulation are large and thus stable time averages are only attained by sampling many rotations. The kinetic energy in the differential rotation and in the convective motions are two orders of magnitude higher than that in the meridional circulation (Brun & Toomre 2002). As a result, small fluctuations in the convective motions and differential rotation can lead to major variations in the circulation. Some of the helioseismic inferences suggest the presence of single cell circulations, which are at odds with our multi-cell patterns. However these inferences vary from year to year, and there is recent evidence for double-cell structure in the circulations observable in the near-surface shear layer, but only in the northern hemisphere as the current solar cycle advances (Haber et al. 2002). From a careful analysis of the angular momentum transport in our shell we have deduced that the slow pole behavior seen in case $`E`$ seems to come about from a relatively weak meridional circulation at high latitudes. This permits a more efficient extraction of angular momentum by the Reynolds stresses from the high latitudes toward the equator in yielding the interesting differential rotation profile that is achieved.

## 4 Incorporating some Magnetic Fields

We now turn to consider the influence of magnetic fields both upon the convection in our deep shell and upon the angular velocity profiles that can be maintained. Early attempts to explain the 22-year solar cycle considered the possibility that the solar dynamo operated within the bulk of the convective envelope (Gilman 1983, Glatzmaier 1985), but such approaches failed because strong magnetic fields could not be stored efficiently within the unstable stratification of the convection zone. The dynamo periods were too short (of the order of 1 year) and the poloidal fields were found to propagate poleward, at variance with observations. More recently, Parker (1993) has proposed an interfacial dynamo model seated in the stable tachocline. Magnetic field is still generated within the bulk of the convection zone ($`\alpha `$-effect), but is pumped downward into the stable layer via overshooting turbulent plumes, to be stretched there into large-scale toroidal structures by the strong shear of the tachocline ($`\omega `$-effect). When the amplification of the toroidal magnetic field is great enough, the structures (called magnetic flux tubes in mean field models) become magnetically buoyant and rise upward through the convective envelope. The strongest of those structures emerge in the photosphere as bipolar magnetic arcades, whereas the weaker ones are recycled within the convective zone. This leads to the crucial natural cycle of poloidal to toroidal exchange, i.e. $`B_{pol}B_{tor}B_{pol}`$. However, many aspects of these essential ‘dynamo building blocks’ remain to be demonstrated through nonlinear 3–D calculations. At present it is not feasible to simulate self consistently all the processes operating together, and thus one needs to concentrate on individual components (magnetic generation, pumping, shearing, buoyant rising). One important ingredient in the interfacial dynamo scenario is the ability of the convective motions to generate and sustain magnetic fields in the bulk of the zone (Cattaneo & Hughes 2001). We thus have evaluated some conditions in 3–D compressible convective shells for which such a dynamo threshold can be realized. We further wished to identify the maximum nonlinear amplification of the magnetic field that can be sustained by the convective motions without destroying the strong angular velocity contrasts previously attained.

### 4.1 Determining the Dynamo Threshold in Convective Shells

We have conducted three MHD simulations (named $`M1`$, $`M2`$ and $`M3`$) started from a solution not as turbulent as case $`E`$ but exhibiting a similar angular velocity profile $`\mathrm{\Omega }`$. We then introduced a small seed dipolar magnetic field and let the simulations proceed. Figure 3 shows the magnetic energy evolution for three values of the magnetic diffusivity $`\eta `$ (i.e. 2, 1.6 & 1 $`\times 10^{12}`$ cm2/s). We note that over more than 3000 days (corresponding to several ohmic decay times) the two lowest diffusive cases $`M2`$ and $`M3`$ exhibit a sustained magnetic energy (ME), the levels of which depend on $`\eta `$. The other case $`M1`$ is clearly decaying, since the rate of generation of magnetic fields could not compensate for the rate of destruction by ohmic diffusion. The dynamo threshold seems to be around a magnetic Reynolds number $`R_{\mathrm{em}}=v_{\mathrm{rms}}D/\eta `$ of $`300`$ or $`\eta 1.6\times 10^{12}`$. This is about 25% higher than in a progenitor incompressible simulation (Gilman 1983).

Further, we find that the kinetic energy (KE) in model $`M3`$ has been reduced by about 40% compared to its initial value. In this case ME has grown to reach a value of 7% of KE. The increase in ME started to influence the total amount of KE contained in the shell when ME reached roughly 0.5% of KE in after about 600 days of evolution. The early exponential growth of ME in case $`M3`$ corresponds roughly to the first 600 days after which the nonlinear feedback of the Lorentz forces on the flow begins to be felt. For case $`M2`$, ME is still small enough (i.e. $`0.1\%`$) even after 4000 days for the convective motions to only be mildly affected by the Lorentz forces. This is most clearly seen in comparing the kinetic energy time trace for cases $`M2`$ and $`M3`$. The magnetic energy stored in the toroidal field within the bulk of the convective zone is roughly an order of magnitude greater than that of the poloidal field.

Figure 4 displays a 3–D rendering of the radial velocity field and of the radial and longitudinal magnetic fields achieved in our convective shells. The radial component of the magnetic field is found to mainly concentrate in the downflow lanes, having been swept away from the center of the convective cells by the broad upflows. Both polarities coexist at the downflow network interstices. The radial component $`B_r`$ possesses finer and more intricate structures than either $`V_r`$ or $`B_\varphi `$, exhibiting many swirls. Substantial magnetic helicity is present, involving complex winding up of the toroidal magnetic fields along their length, with both polarities interchanging their positions. The toroidal magnetic fields also possess the greater spatial scales, having been stretched by the gradients in angular velocity especially in the equatorial region. Some features could resemble magnetic flux tubes, although they are short lived and less concentrated.

### 4.2 Effect of the Lorentz Force

With fairly strong magnetic fields sustained within the bulk of the convection zone in case $`M3`$, it is to be expected that the differential rotation $`\mathrm{\Omega }`$ will respond to the feedback from the Lorentz forces. Figure 5 shows the time averaged angular velocity achieved in case $`M3`$. The main effect of the Lorentz forces is to extract energy from the kinetic energy stored in the differential rotation. The reduction of KE contained in the angular velocity is of the same order as the decrease seen in the total KE, i.e. 40%. As a consequence, the angular velocity contrast $`\mathrm{\Delta }\mathrm{\Omega }`$ from 60 to the equator drops by $`30\%`$ in case $`M3`$, going from 140 nHz (or 34% compared to the reference frame $`\mathrm{\Omega }_o`$) in the hydrodynamic case to 100 nHz here (or 24%). Nevertheless, the angular velocity in case $`M3`$ remains in reasonable agreement with the solar profile. The reduction of the latitudinal contrast of $`\mathrm{\Omega }`$ can be attributed to the poleward transport of angular momentum by the Maxwell stresses (the mean magnetic fields having a negligible contribution). Now the Reynolds stresses again need to balance the angular momentum transport by the meridional circulation, the viscous diffusion and the Maxwell stresses. This leads to a less efficient speeding up of the equatorial regions. Since ME is only 7% of KE in case $`M3`$, the Maxwell stresses are not yet the main players in redistributing the angular momentum. We have found that a value of ME above about 20% of KE leads to a significant magnetic braking effect on the differential rotation. Had the simulation been restarted with a stronger initial magnetic field $`B_0`$, $`\mathrm{\Delta }\mathrm{\Omega }`$ could drop by 90% in less than a few hundreds days, thus being at variance with helioseismic findings. Thus the feed back from the Lorentz forces can give an upper limit on the maximum sustainable magnetic field for a given rotation profile. In case $`M3`$ it seems that the limit is a few thousand gauss but we have to be cautious since this computation does not allow for expulsion of magnetic helicity or pumping of magnetic field into the radiative zone.

### 4.3 Magnetic Activity and Possible Cycles

I now wish to discuss possible differences between active stars and ones with solar-type cycles of magnetic activity. The interfacial dynamo invoked in this work to explain the solar cycle require a strong rotational shear located in the stable tachocline. Without such shear, the star is more likely to undergo local magnetic dynamo action directly associated with helical MHD turbulent convective flows (i.e the $`\alpha `$-effect). Yet, any convection zone under the influence of rotation is likely to have a more or less pronounced differential rotation (cf. Brun, Browning & Toomre 2003 for core convection and dynamo processes in A-type stars). Thus some $`\omega `$-effect can also be expected within a convection zone, except that the strongest magnetic fields present in the unstable convective zone will certainly be smaller than the ones amplified in the stable tachocline (again because of buoyancy considerations). Thus, it is unlikely that purely convective stars such as M dwarfs or pre-main sequence stars will exhibit cycles, but instead may display irregular activity (the strength of which will likely depends on the rotation rate or Rossby number $`R_o`$).

## 5 Conclusions

Our simulations of convection in 3-D spherical shells are describing how turbulent convection under the influence of rotation can achieve a well organized differential rotation. Transport of angular momentum by tilted turbulent plumes generate Reynolds stresses leading to an equatorial acceleration, such as seen in the Sun. When we introduce a seed magnetic field and let the simulation evolve until equilibration we find solutions where it is possible to sustain both a magnetic field and a solar-like differential rotation. The amount of magnetism that can be sustained in a convective shells for a given rotation profile gives us some clues on how solar magnetism is working. Strong energetic event like CME’s and the pumping of magnetic energy in the stable radiative zone certainly plays a significant role in regulating the intricated and coupled interplay of convection, rotation and magnetism. From a more general stellar point of view, it is clear that convection under the influence of rotation and magnetism will lead to significant activity but not necessarily to cycles. For the latter, a stable shear region such as the solar tachocline seems to be a crucial element in order to close the loop of magnetic field generation and organization.

I thank Mark Miesch, Annick Pouquet, Juri Toomre and Jean-Paul Zahn for useful discussions and collaborations and Jean Arnaud for its invitation to the workshop. This work was partly supported by NASA through SEC Theory Program grant NAG5-8133 and by NSF through grant ATM-9731676. The simulations with ASH were carried out with NSF PACI support of the San Diego Supercomputer Center (SDSC). Much of the analysis of the extensive data sets has been done in the Laboratory for Computational Dynamics (LCD) within JILA.