text:
# modeling study on the validity of a possibly simplified representation of proteins

## I INTRODUCTION

Protein folding is a well known complicated and highly cooperative dynamic process due to the heterogeneity in proteins (e.g., see Refs. and references therein). Much effort has been made by considering minimalist models with a few types of amino acid residues to simplify the natural set of residues (of 20 types) for better physical understanding  and practical purposes . In these models the compositions are much simpler than the real ones. The simplest reduction is the so-called hydrophobic and polar ones with only two groups of residues (each group including some types of residues) by considering the main driving force, the hydrophobicity . Furthermore, these two groups actually are simplified as two effective monomers or letters, namely H and P, which is known as the HP model . The studies of such a model enable people to understand some fundamental physics and mechanism of protein folding. However, as argued in a number of studies (see  and reference therein), the HP model may be too simple and lacks enough consideration on the heterogeneity and the complexity of the natural set of residues, such as the interactions between the residues . Moreover, the minimal sets of residues for protein design suggested by biochemical experiments  seem unfavorable to those with only two types of residues since a small number of types obviously introduces the homopolymeric degeneracy. What is the suitable simplification for natural proteins, or how many types of residues are necessary for reproducing some useful structures and for a simplified representation of protein sequence characteristics? These are not well understood.

Recently, Riddle and et al. made an exciting approach to the problems mentioned above experimentally . By using combinatorial chemistry along with a screening strategy, they searched and found out a subset of the natural amino acids that can be used to construct a well-ordered protein-like molecule consisting of $`\beta `$-sheets. This subset contains five amino acids , namely, isoleucine, alanine, glycine, glutamic acid and lysine, which are simply represented as I, A, G, E and K. Although about $`30\%`$ of the residue sites in some of their sequences are not encoded by the five-letter palette, rather than by nine other types of residues, the sequence complexity of the protein is largely simplified. These uncoded sites are due to their direct involvement in binding a proline-rich signal peptide. \[if the binding was compromised, the protein would not show up in the screening array .\] They also found that subtractions of the five-letter code to a three-letter one may destroy the structural rebuilding. As argued by Wolynes, this experiment shows the possibility of simplification for the protein sequence complexity, and that some complexity with five-letter code, but not three-letter code, might be still needed based on the landscape ideas. This experiment extends the search for minimal solution of the simplified representation of complex protein architecture. However, it is not clear whether the suggested 5-letter code is valid in general and feasible for elucidating characteristics of real proteins with 20 kinds of amino acids. That is, is the 5-letter code based on a specific protein generally workable for other natural proteins?

In this work, we address the questions mentioned above at the level of a lattice model with some contact-potential-form interactions. Based on the statistical and the kinetic characteristics of the folding, and on the thermodynamic stability of the ground states of some reduced sequences, we study the validity of two reductions, namely one with the five-letter palette I, A, G, E and K  and the other composed by two types of amino acids as an example of HP-like model. We find that the five-letter reduced sequences display similar native structural features and folding kinetic behavior to the optimal 20-letter ones. Differently, for the two-letter case, our results show its deficiency to act as a good reduction for general representation of proteins. The five-letter reduction may be a suitable description for simplifying the sequences complexity of 20-letter model chains which may relate to the natural proteins. A detailed analysis on correlation between the sequences with 20-letters and their five-letter substitutes shows the source of such validity of the five-letter code.

## II MODELS

In general, a successful reduction implies that an original sequence with 20 types of residues can be represented by a new sequence with fewer types. The reduced sequence should have the main statistical and structural characteristics, as well as the kinetic accessibility, as that of the original one, such as the basic residue components of the original sequence, the ground state or the native conformation, etc. That is, after the reduction, not only should the statistical characteristics of the energy spectrum (the existence of a large ground state gap) be maintained, but also the folded structure should be kept the same as that of the original one (see Fig.1). Such a reduction is regarded as a good one. Nevertheless, if the energy gap is diminished, or the ground state of the reduced sequence is degenerate (or deviates from that of the original one), the thermodynamic stability and the kinetic characteristics are altered. Consequently, the reduction is believed to be a bad one. Therefore, to explore the validity of simplified reduction means to verify it being a good one or not.

Let us start by considering a model protein which is assumed to be a chain with connected residues in a cubic lattice (Fig.2a). A self-avoided arrangement of the model chain on the lattice is generally noted as a structure or conformation of the model protein. The conformation of the model chain is characterized by the set $`\{𝐫_𝐢\}`$, the spatial position of the $`i`$-th residue in the chain, and its sequence by $`\{s_i\}`$, the set of various types of residues assigned along the chain. An example of such a chain is shown in Fig.2(a), with $`s_1=V`$, $`s_2=P`$, and $`s_3=V`$, and so on, and the chain length is $`L=27`$, i.e., $`27`$ residues. The energy of the chain is determined by its sequence and conformation,

$$E(\{s_i\},\{𝐫_𝐢\})=\underset{i>j}{\overset{L}{}}B(s_i,s_j)\mathrm{\Delta }(𝐫_𝐢,𝐫_𝐣),$$
(1)
where $`B(s_i,s_j)`$ is the contact potential between the $`i`$-th residue (with residue type $`s_i`$) and the $`j`$-th residue (with residue type $`s_j`$), and the function $`\mathrm{\Delta }(𝐫_𝐢,𝐫_𝐣)=1`$ when the $`i`$-th and the $`j`$-th residues are spatial, not sequential, neighbors and $`\mathrm{\Delta }(𝐫_𝐢,𝐫_𝐣)=0`$ otherwise. To study the reduction, a number of such sequences, say $`N=100`$, made of the natural set of amino acids with 20 letters are all well optimized to a same native structure \[the one shown in Fig.2(a)\] using the methods proposed by Shakhnovich et al. All these sequences are referred as “original” ones. These original sequences then are reduced according to a pre-selected five-letter reduction scheme as follows : group-I (with residues C, M, F, I, L, V, W and Y), group-II (A, T and H), group-III (G and P), group-IV (D and E) and group-V (S, N, Q, R and K), with a representative residue for each group as I, A, G, E and K, respectively. That is, an original sequence is reduced into a sequence of five-letters by substituting (or replacing) each letter with its representative letter. In such a reduction, 20 types of residues are divided into five groups according to their interaction characteristics (following the Miyazawa and Jernigan (MJ) matrix). Each group contains some residues which interact with others in a similar way . For example, group-I includes all the hydrophobic residues, and the other four contain the residues with polar features. Obviously, this five-group simplification takes into account of more heterogeneities for the protein, and also considers more detailed differences between the polar residues than the two-letter HP model . Then we take Baker’s 5-letter alphabet, I, A, G, E and K, as the best representative letters for five groups based on a physical reason, but not an arbitrary choice (A more detailed description and discussion see Ref.).

The procedure to reduce the original sequences to two-letter cases is similar as above. Here, we consider group $`H`$ including residues C, M, F, I, L, V, W and Y, and group $`P`$ including the rest, and we pick the residues I and A as the representatives of group H and P, respectively. Some other choices, such as I and E, are also checked. It is found that these choices give out basically the same results.

## III RESULTS AND DISCUSSIONS

Does the reduced sequence with five letters have the same folded native structure as that of the original one shown in Fig.2(a)? Let us make a statistics on the ratio of successful reduction ($`RSR`$), for the reduced sequences. Because of the main attractive feature of the MJ interaction matrix, the native state generally has a structural motif with the form of maximally compact structure. Besides, for the 27-mer model chain we studied, the compact conformations are generally understood as the $`3\times 3\times 3`$ cube-shape structures. Thus, the statistics on the RSR can be found out by enumerating all possible compact structures in the case of a $`3\times 3\times 3`$ lattice . Interestingly, almost $`86`$ reduced sequences (out of $`100`$) with five-letters are “foldable” since they have unique ground states. Among the foldable sequences, $`74`$ keep the same native structure of the original sequences, termed as “folded”, and $`12`$ have different folded states, termed as “ground states deviation” \[see Fig.2(b)\]. This means that this five-letter strategy basically reproduces the original native structure with $`RSR`$ of $`74\%`$. In addition, there are $`14`$ sequences that do not have unique ground states, termed “unfoldable”.

As a comparison, the $`RSR`$ is also evaluated for the HP reduction. It is found that there are a few foldable sequences, namely, three reduced sequences (out of $`100`$) but with “ground states deviation”. That is, these three sequences fold to unique ground states which are different from that of the native state of the original sequences, i.e., the $`RSR`$ is $`0`$. However, the existence of the unique ground states basically coincides with the result in Ref. where there is only 4.57% of the sequences having unique ground states. Physically, it results from the homopolymeric degeneracy as argued in Ref. and found in other lattice simulations. The biochemical experiments also approve this deficiency, since it is found that two types of residues are far insufficient for the rebuilding of four-helix bundle proteins in protein engineering experiments and in de novo designs .

Differently, the five-letter case is quite good in overcoming the homopolymeric tendency. The similarity between the five-letter case and the 20-letter case with a large value of $`RSR`$ indicates the validity of the five-letter strategy in simplifying the natural proteins. It is noted that the selection of the optimal original sequences and their native structures is rather arbitrary without any bias on interactions and compositions. Several studies on some other sequences with different native structures present similar results, which further implies that this five-letter reduction scheme is universal for different structures. Therefore, the simplified depiction of the Riddle et al.’s five-letter substitutes may catch basically the physics of natural proteins in their structures.

What kinetic characteristics then do the five-letter sequences show? These are further studied by the standard lattice Monte Carlo simulations. The thermodynamic quantities, such as heat capacity $`C_v`$ of the model chains, are calculated based on the histogram algorithm over a collection of samples of conformations. In our simulations, the native structures found in the above enumerations over various compact structures are checked. It is found that the checked structures have uniquely minimal energies and frequently appear in the sampling for these structures. These ensure the thermodynamic stability and kinetic accessibility. It approves their characteristics as native structures found in our enumerations. Therefore, the statistics on the $`RSR`$ based on these kinetic Monte Carlo simulations shows basically the same value five-letter reduction as that obtained by enumerations. Similar results have also been obtained for the HP reduction.

Furthermore, we use the foldability $`\sigma =|T_\theta T_f|/T_\theta `$ proposed by Klimov and Thirumalai as an indicator of the kinetic characteristic for the system. Here $`T_\theta `$ indicates the temperature corresponding to the peak of heat capacity of a model chain, and $`T_f`$ marks the structural transition to the well-ordered native structure during the folding processes, which corresponds to the maximum of the fluctuation of structural overlap function $`\chi (T)`$. A small $`\sigma `$ is generally related to a fast folding with few kinetic traps. On the other hand, a large $`\sigma `$ ($`1`$) may result in many competitive local minima, and in this case a protein can not fold in a biologically relevant time scale. Thus, the factor $`\sigma `$ is argued as a criterion of the foldability although there is a controversy on its definition. The correlation between the factor $`\sigma `$ and folding ability has been illustrated in many cases (see Ref. and references therein). Here, we use the factor $`\sigma `$ to characterize the foldability of our model chains. Note that we do not monitor the actual folding time of the model chains in our simulations, and our following discussions on kinetics depend on the presumed relationship between the kinetics and thermodynamics, i.e., the foldability $`\sigma `$ proposed in Ref.. A relevant discussion on such a relationship is recently made by Dinner et al .

We then calculate the foldability $`\sigma `$ for the original sequences and their five-letter and two-letter substitutes, respectively (see Fig.3, the solid circles). One can see that for the five-letter substitutes the foldability is $`\sigma =0.46`$ which is smaller than $`0.6`$, a critical value of slow folding sequences . This suggests that the native state is accessible kinetically for the five-letter reduction, and shows the similar dynamics of two kinds of sequences, namely, the original ones and the reduced ones. Differently, for the two-letter case the foldability is $`\sigma 1`$ much larger than 0.6, which means little kinetic accessibility of its native structure. In addition, some independently optimized sequences (i.e., some 20-letter sequences optimized with the five-letter or two-letter alphabet, respectively) are also studied. The values of $`\sigma `$ are found to be 0.25 and 0.75 for the five-letter and two-letter sequences, respectively, which shows clearly the improvement of the kinetic accessibility  due to the optimization for each sequence. Moreover, a similar exponential decreasing tendency can be seen (see Fig.3, the open diamonds). As a result, both exponential tendencies of decrease clearly imply the validity of a five-letter strategy.

It is worth noting that in our simulations, a large number of extended conformations are explored, not only the most compact conformations within the $`3\times 3\times 3`$ cubic lattice. These kinetic folding simulations may be important for the case of longer chains since the non-compact conformations have more effects on the folding features of the model proteins . Fortunately, the competition of non-compact structure as a native candidate may not be so serious in our case (see the remarks at the end of this sections). Certainly, for a more realistic and detailed study, we may need to make more extensively kinetic folding simulations on longer chains, which deserves further work. However, we may believe that the basic physics is the same because many main characteristics of the folding have been found from the lattice model of proteins previously .

Now, let us discuss the correlation between the original sequences with 20-letter code and their substitutions with five-letter code. As we all know, there is $`m^L`$ possible sequences for the proteins with $`m`$ types of residues and length $`L`$. Can the optimized part of $`20^L`$ sequences be successfully mapped to $`5^L`$ substitutes as mentioned-above? To answer this question, we make a detailed analyses on the sequence set $`S_{20}^0=\{S_{20}^1,S_{20}^2,S_{20}^3,\mathrm{}\}`$ in which all their substitutes with the five-letter reduction are the same, i.e., a single sequence $`S_5^0`$. We randomly select a number of compact structures, say 50, from the five-letter spectrum as targets for design, i.e., optimizing the sequences $`S_{20}^0`$ to lower energy .

In practice, from a sequence $`S_5^0`$, we produce a number of sequences $`S_{20}^1`$, $`S_{20}^2`$, $`\mathrm{}`$, $`S_{20}^{30}`$, randomly following the reduction scheme. That is, each of I,A,G,E, and K is randomly substituted by a residue belonging to its corresponding group, e.g., A by A, T, or H randomly. Then for a certain target structure, we process the design for such a sequence $`S_{20}^i`$ by randomly exchanging their positions of two residues which belong to the same group (keeping the five-letter reduced sequence $`S_5^0`$ unchanged). The new sequence is accepted when the energy increment $`\mathrm{\Delta }E<0`$, otherwise accepted with a probability $`P=exp(\mathrm{\Delta }E/T)`$ with $`T=0.1`$. Until $`10^7`$ new sequences are reached, the optimized sequence is found out for which its energy is the lowest one. The target is designable if its energy is the lowest one, otherwise undesignable.

We find three cases: 1) Undesignable compact structures (UDS): Some compact structures never do become the native states of the set $`S_{20}^0`$ even after very long time optimization and very slow annealing (see Fig.4). 2) Designable compact structures with small energy gaps for their ground states (DSSG): some compact structures are designable and the related sequences fold exactly into these structures. However, the energy gaps above their ground states are quite small, generally less than 0.5 (see Fig.4). According to the features of the protein energy spectrum, such small energy gaps result in an unstable ground state , and also many traps in the energy landscape . Thus these proteins may fold rather slow and may be unstable in their folded states. 3) Designable compact structure with large energy gaps for their ground states (DSLG): there exists a specific compact structure, namely, the native structure of the sequence $`S_5^0`$. Taking it as a design target, the designed sequences, from all of $`S_{20}^0`$, not only fold to it, but also have large energy gaps and show good foldability (see Fig.4). This implies that a set of optimal sequences $`S_{20}^{}`$ that have the same five-letter reduced sequence shows stable folding and is included in the set $`S_{20}^0`$, $`S_{20}^{}S_{20}^0`$. Thus, if two optimized sequences $`S_{20}^i`$ and $`S_{20}^j`$ have the same 5-letter reduced sequence $`S_5^0`$, they will behave the same in folding behavior as that of $`S_5^0`$. \[If sequences with 20 letters are optimized without the condition of keeping $`S_5^0`$ unchanged (see the preceding discussions), most of the optimized sequences (about 74%) will have the same folding behavior after the substitution (see Fig.2(b))\]. In other words, the complexity of the sequences $`S_{20}^0`$ can be simplified via the reduction since both the $`S_5^0`$ sequence and $`S_{20}^0`$ sequences behave almost the same in the aspects of folding kinetics and folded structure. Thus by the reduction the sequence space (with 20 letters) can be reduced to many sets, and each set has a single five-letter sequence and a favorable native structure. In this level, the correlation of multiple-to-one between a set of 20-letter sequences and a five-letter sequence is well established and ensures the similarity of two sequence spaces. Finally, it is noted that other five-letter reductions do not present such simplification from the original sequences, and have a small (large) value of $`RSR`$ ($`\sigma `$) since they are not the best representative letters for groups .

Finally, let us remark on our modeling study. It is argued that the enumerations over only the compact structures are questionable for searching the native state of a lattice chain. Indeed, this is true for the usual HP model (with $`E_{HH}=1`$ and $`E_{HP}=E_{PP}=0`$)  and that used by Li, et al. (with $`E_{HH}=2.3`$, $`E_{HP}=1`$, and $`E_{PP}=0`$) . In those models, there are some binding interactions between residues, such as $`E_{PP}=0`$, that are as weak as that between a residue and solvent molecules. This may mean that there is little difference in energies between the compact structures and the non-compact ones in those models , which introduces some limits in enumeration analysis. However, the interactions for our chains, not only for 20-letter chains, but also for the five-letter cases and two-letter cases, are taken from the MJ matrix directly, which has an obvious attraction between any pairs of residues. Thus, there is a large energy penalty for the extended structures with fewer bonds, i.e., an energy bias towards maximally compact structure. It improves the stability of the compact structures in our analysis. This makes our search for the native structures by enumerations over the set of compact structures feasible. In addition, the lattice chains with 27 monomers are shorter than those analyzed in Ref.. As a result, we believe that in our cases the determination of the native structures based on the enumerations over compact structures is appropriate, though there is a mixing between the energy spectrum (the energy levels above the native state) of the compact structures and those of the non-compact structures, which may affect the kinetic properties. In practice, the native structures obtained by enumerations in our work are verified by the Monte Carlo kinetic simulations. That is, for our cases, especially the five-letter chains, we have not found any other states with lower energy than the native structures during very-long-time simulations. Besides, using enumeration over compact structures is more effective than the kinetic simulations, which enables us to make an extensive search over the sequence space on the structural rebuilding problem, as interested in our work.

## IV SUMMARY

In this work, the five-letter code for simplification of proteins suggested by Riddle et al. is found to be valid not only keeping the features of the energy spectrum but also reproducing the native structure of the sequences with the natural set of residues. The kinetic resemblance between the sequences composed with 20 letters and five-letters also implies a similarity of their funnel characteristic. This illustrates that the lattice model of proteins with natural set of residues can be re-constructed with a smaller set of residues, maybe five types of residues. Thus our study suggests that the five-letter code can act as a suitable simplification considering more heterogeneity of the model proteins, and may encode the main information for model proteins. In addition, since that the folding kinetics in real proteins and lattice models is believed to be similar in many aspects , our simulations on the lattice chains may provide some understanding on the simplification of real proteins. As a comparison, the unfavorable reductions with two types of residues may indicate that the proteins with a natural set of residues are too complex to be simplified by the two-letter representations. That is, the minimal solution of protein simplification could not be as simple as the HP models, at least in lattice level. This actually coincides with the conclusion in several protein design experiments . Finally, we note that an exact and detailed grouping of residues for the simplification representation of proteins still needs more theoretical and experimental explorations.

## V ACKNOWLEDGMENTS

* W.W. acknowledges support from the Outstanding Young Research Foundation of the National Natural Science Foundation of China (Grant No. 19625409). We thank C. Tang for a critical reading of the manuscript. The numerical simulations in this work were done on the SGI Origin 2000 in the National Laboratory of Solid State Microstructure, Nanjing University.


format: text

text:
# Untitled Document

Thermal expansion and the equation of state of Ir and Rh

By M. I. Katsnelson, M. Sigalas, A. V. Trefilov and K. Yu. Khromov

Russian Research Centre “Kurchatov Institute”, Moscow 123182, Russia
Institute of Metal Physics, Ekaterinburg 620219 ,Russia
$``$Ames Laboratory, Iowa State University, Ames, IA, 50011, USA

Abstract

The simplest anharmonic characteristics of Ir and Rh are discussed in the framework of a previously developed simple pseudopotential model which describes the elastic moduli, phonon spectra and the lattice heat capacity in the harmonic approximation of these metals succesfully. The microscopic Gruneisen parameters, the dependences of the elastic moduli on pressure, the coefficient of thermal expansion and the equations of state at the finite temperatures have been calculated. The ab initio  calculations of the energy-band structure and the equation of state for Ir at $`T=0`$ have been done to test the model for adequacy at high pressures. The values of different contributions (zero-point oscillations, quasiharmonic, etc.) in the considered thermodynamic characteristics of Ir and Rh are discussed.

§1. Introduction

As it is stated in the literature (Gornostyrev et al 1994, Ivanov et al 1994, Katsnelson et al 1996, hereinafter referred to as I), Ir and its analogue Rh are set off from other FCC metals because of some distinctive features of their physical properties: unusual deformation-induced failure, peculiar defect structure, specific temperature dependence of the effective Debye temperature, etc. Moreover, the situation is not trivial because in the case of Ir, unlike the most, if not all of the other FCC metals a rather broad range of properties (elastic moduli, phonon spectra and lattice heat capacity) can be described in terms of a simple pseudopotential model (see I and references therein). Even in $`sp`$-metals the situation is less favorable: to describe the elastic moduli of Ca allowance must be made for the singular contributions to energy which result from the proximity of the Fermi surface and the faces of the Brillouin zone (Katsnelson et al 1990), whereas for Al it is necessary to take into account the contributions of three body forces. Nevertheless, in describing the phonon spectra of Al it is impossible to attain an accuracy identical to that for Ir (Browman et al 1971, Fomichev et al 1990). In Ir the contribution of three body interactions proves to be small (Greenberg et al 1990). As a whole, it turns out that surprising as it may seem, Ir is probably exceeded in the accuracy of the psedopotential description of the lattice properties only by alkali metals. This gives hope that the theory can be used to predict those lattice properties of Ir for which the experimental data are scanty or absent, such as thermal expansion or the equation of state at different temperatures. To be sure, the pseudopotential model is oversimplified. Therefore, where it is possible, the model predictions should be tested against either the experiment or the results of the ab initio  band calculations. Despite the fact that no simple model can completely substitute the ab initio  calculations, to obtain a sufficiently reliable physical information about the complex properties of metals even with the help of a model approach appears to be of interest. At the same time, fully ab initio  calculations of the phonon contributions to the thermodynamical properties of metals seem for the present to be impossible (or, in any case, are rather cumbersome). Therefore, the use of the pseudopotential model for these purposes, where it is possible seems to be justified. This paper presents the results of the calculations of the coefficient of thermal expansion at finite temperatures, and the dependence of the elastic moduli on pressure in terms of the pseudopotential model developed previously for Ir and Rh (see I). To understand to what extent the simple model used is reliable at high pressures the electron energy spectrum of Ir and its equation of state at $`T=0`$ were calculated ab initio . The results of the calculations show that the equation of state calculated in the framework of the pseudopotential model seems to be sufficiently reliable.

§2. Thermal expansion and Gruneisen parameters of Ir and Rh

In calculating the thermal expansion and the Gruneisen parmeters $`\gamma (T)`$ at finite temperatures as well as the dependencies of the elastic moduli on pressure $`B_{ik}(p)`$ use was made of the pseudopential model constructed in I and giving the optimal description of the phonon spectra. The results of the calculations in terms of an earlier model (Ivanov et al 1994) are given for comparison. The local pseudopotentials $`V_{ps}(q)`$ used in the two models are closely similar. However a slight change in the form of $`V_{ps}(q)`$ was made in I to describe better the phonon spectra. The screening function was taken in both cases in the Geldart and Taylor (1970) approximation (see, in more detail, Greenberg et al 1990). The temperature dependence of the volume per atom $`\mathrm{\Omega }(T)\mathrm{\Omega }_0+\mathrm{\Delta }\mathrm{\Omega }(T)`$ ($`\mathrm{\Omega }_0=\mathrm{\Omega }(T=0)`$), the coefficient of thermal expansion

$$\alpha _p=\frac{1}{\mathrm{\Omega }}\left(\frac{\mathrm{\Omega }}{T}\right)_p$$
$`(1)`$
and the Gruneisen parameter $`\gamma (T)`$ are determined by the following expressions (Vaks et al 1978):

$$\frac{\mathrm{\Delta }\mathrm{\Omega }(T)}{\mathrm{\Omega }_0}=\frac{p_{ph}(T)}{B_T},$$
$`(2)`$
$$\alpha _p=\frac{1}{\mathrm{\Omega }B_T}\underset{\lambda }{}\left(\frac{\mathrm{}\omega _\lambda }{T}\right)^2N_\lambda (1+N_\lambda )\gamma _\lambda $$
$`(3)`$
$$\gamma (T)=\frac{\mathrm{\Omega }B_T\alpha _p}{C_V(T)}=\frac{{\displaystyle \underset{\lambda }{}}\left({\displaystyle \frac{\mathrm{}\omega _\lambda }{T}}\right)^2N_\lambda (1+N_\lambda )\gamma _\lambda }{{\displaystyle \underset{\lambda }{}}\left({\displaystyle \frac{\mathrm{}\omega _\lambda }{T}}\right)^2N_\lambda (1+N_\lambda )}$$
$`(4)`$
where

$$p_{ph}=\frac{1}{\mathrm{\Omega }}\underset{\lambda }{}\mathrm{}\omega _\lambda N_\lambda \gamma _\lambda $$
$`(5)`$
is the phonon pressure; $`\lambda 𝐪\nu `$, q is the phonon wave vector (sweeping the Brillouin zone); $`\nu `$ is the branch number; $`\omega _\nu `$ is the phonon frequency;

$$\gamma _\lambda =\frac{\mathrm{ln}\omega _\lambda }{\mathrm{ln}\mathrm{\Omega }}$$
$`(6)`$
are the microscopic Gruneisen parameters;

$$N_\lambda =\frac{1}{\mathrm{exp}{\displaystyle \frac{\mathrm{}\omega _\lambda }{T}}1},$$
$`(7)`$
$`B_T`$ is the bulk modulus at a constant temperature; $`C_V(T)`$ is the phonon heat capacity at a constant volume. The microscopic Gruneisen parameters are calculated by the formula

$$\gamma _\nu (𝐪)=\frac{\mathrm{\Omega }}{2\omega _{\nu 𝐪}^2}\underset{\alpha \beta }{}(e_{\nu 𝐪}^\alpha )^{}\frac{D_{\alpha \beta }(𝐪)}{\mathrm{\Omega }}e_{\nu 𝐪}^\beta $$
$`(8)`$
where $`D`$ is the dynamic matrix; $`𝐞_{\nu 𝐪}`$ are the polarization vectors; $`\alpha `$ and $`\beta `$ are the Cartesian indices. In our calculations allowance was also made for the quasiharmonic corrections:

$$\alpha =\alpha _1\left[1\left(\frac{\mathrm{\Delta }\mathrm{\Omega }}{\mathrm{\Omega }}\right)_1\frac{\mathrm{ln}B}{\mathrm{ln}\mathrm{\Omega }}\right]$$

$$\frac{\mathrm{\Delta }\mathrm{\Omega }}{\mathrm{\Omega }}=\left(\frac{\mathrm{\Delta }\mathrm{\Omega }}{\mathrm{\Omega }}\right)_1\left[1+\left(\frac{\mathrm{\Delta }\mathrm{\Omega }}{\mathrm{\Omega }}\right)_1\left(1\frac{\mathrm{ln}B}{\mathrm{ln}\mathrm{\Omega }}\right)\right]$$
$`(9)`$
where index “1” means the corresponding expression taken in the lowest order of anharmonicities. The procedure of calculating the thermal expansion of metals and its related quantities was discussed in more detail in Vaks et al 1978.

The results of the calculations are presented in figs. 1–4 and table 1. First of all, it should be noted once more that $`\gamma _\nu (𝐪)`$ vary essentially over the Brillouin zone. This indicates again that the Gruneisen approximation $`\gamma _\nu (𝐪)=const`$ is completely unsuited for actual systems. As the temperature increases the macroscopic Gruneisen parameter $`\gamma (T)`$ reaches quickly (at $`T0.1\mathrm{\Theta }`$, where $`\mathrm{\Theta }`$ is the Debye temperature) a constant high temperature value. Alkali metals exibit the analogous behavior (Vaks et al 1978). However, unlike alkali metals having $`\gamma 1`$, in Ir and Rh, as in majority of transition metals, $`\gamma 2`$. Note also that the calculated values of $`\gamma (T)`$ are rather highly sensitive to the form of the pseudopotential used (compare the different curves in figs 3b and 4b).

§3. Equation of states. Used formulas and approximations

The equation of state for a metal, that is, the dependence of the pressure $`p`$ on $`\mathrm{\Omega }`$ and $`T`$ is determined by the expression

$$p(\mathrm{\Omega },T)=p_0(\mathrm{\Omega })+p_{zp}(\mathrm{\Omega })+p_{ph}(\mathrm{\Omega },T)+\mathrm{\Delta }p_e(\mathrm{\Omega },T)$$
$`(10)`$
where

$$p_0(\mathrm{\Omega })=\frac{E(\mathrm{\Omega })}{\mathrm{\Omega }},$$
$`(11)`$
$`E(\mathrm{\Omega })`$ is the total energy of the crystal at $`T=0`$;

$$p_{zp}=\frac{}{\mathrm{\Omega }}\underset{\lambda }{}\left(\frac{\mathrm{}\omega _\lambda }{2}\right)=\frac{1}{2\mathrm{\Omega }}\underset{\lambda }{}\mathrm{}\omega _\lambda \gamma _\lambda $$
$`(12)`$
is the pressure relating to the zero-point energy; $`p_{ph}`$ is the phonon pressure (5);

$$\mathrm{\Delta }p_e=\frac{\pi ^2}{6}T^2\frac{N(E_F)}{\mathrm{\Omega }}$$
$`(13)`$
is the temperature dependent contribution to the electron pressure; $`N(E_F)`$ is the density of states at the Fermi level $`E_F`$. We shall neglect the anharmonic contributions having the same order of magnitude as the deviation of the $`\mathrm{\Delta }p_e(T)`$ from the “low-temperature” ($`TE_F`$) expression (13). Generally speaking, the allowance for the contribution of the zero-point vibrations $`p_{zp}(\mathrm{\Omega })`$ to the equation of state in the scheme with a pseudopotential whose parameters are fitted, among other things to the condition $`\mathrm{\Omega }=\mathrm{\Omega }^{exp}`$ at $`p=0`$ requires to refit the parameters. However the results of Vaks et al 1977 show that even in metals as light as Li and Na such refitting is not very essential. In addition, in transition metals, in particular in Ir, the choice of the value for the effective charge is not quite unumbiguous, and strictly speaking, $`Z`$ can be dependent of $`\mathrm{\Omega }`$ ($`Z1`$ in alkali metals). Because we ignore this knowingly more important effect, the results obtained in refitting the pseudopotential parameters with allowance for $`p_{zp}`$ are not presented. The direct calculation shows that that the change in the results with such a refitting is actually small.

In calculating the major contribution to the pressure $`p_0(\mathrm{\Omega })`$ by (11) in the framework of the pseudopotential model for the energy $`E(\mathrm{\Omega })`$, use was made of the expression derived in the second order of the perturbation theory from a pseudopotential with the exchange-correlation energy according to analitical approximation (Perdew and Zunger 1981) of the results of the calculation by the quantum Montne-Carlo method (Ceperly and Alder 1980). For more details see Greenberg et al 1990. The ab initio  calculations of the band structure and the total energy at different $`\mathrm{\Omega }`$ were done by the augmentned plane wave method. The detailed procedure of these calculations is described by Sigalas et al 1992 and Sigalas and Papaconstantopulos 1994.

§4. Equation of state. Results of calculations

The strict calculation of equation of state for a metal by Eq. (10) is rather difficult. Whereas the electron contributions $`p_0`$ and $`\mathrm{\Delta }p_e`$ can be calculated form “the first principles”, the calculation of the phonon contributions $`p_{zp}`$ and $`p_{ph}`$ in the fully “first- principle” approach is a very cumbersome procedure, and no corresponding results are known to us at that time. For this reason, to calculate the phonon contributions to pressure requires to use some model concepts.

We’ll discuss first the results of the $`T=0`$ calculations. Fig. 5 shows the values of $`p_0(\mathrm{\Omega })`$ calculated in the framework of the pseudopotential model and the band approach that we used. It is seen that the difference of the two curves is not too great, at least, at moderate pressures. The reason is that the values of the bulk modulus: 3.76 Mbar calculated with the used variant of the band theory (Sigalas et al 1992) and 2.93 Mbar in the pseudopotential model (see I) do not differ too much.

Table 2 lists the contribution of the zero point vibrations to the equation of state at $`T=0`$. It is seen that this contribution is negligible. The results of the calculation of the phonon pressure in the pseudopotential model are shown in figs. 6 and 7. The temperature dependent electron contribution to the pressure was derived from the results of the band calculations according to (13). It appears to be very small in comparison with the phonon contribution (50 times smaller at the temperature 295 K). The resulting equation of state for Ir at different temperatures is presented in fig. 8. The results for Rh are similar and not presented here. Unfortunately, at present we have no experimental data to be compared with the above results.

§5. Elastic moduli at high pressures

Information about the behavior of the Ir and Rh lattice properties at different pressures can be obtained in the framework of the simple pseudopotential model used here. The extraction of such information from experiments seems today to be difficult. Therefore, we present the results of the corresponding calculations for Ir (fig. 9). The results for Rh are similar.

§6. Conclusion

The reasons why the lattice properties of Ir and Rh are so succesfully described in the framework of the used simple pseudopotential model are not fully understood. However it is believed that if not sufficient, then at least necessary condition for this is the absence of any noticable singularities in the density of electron states in the immediate vicinity of $`E_F`$. The results of the calculations made by us at different pressures (figs. 10,11) show that at not too high pressures this remains the feature of Ir, and, what’s more, the character of the electron structure near $`E_F`$ is unchanged. Therefore, it can be assumed that whatever the reason for the succesful pseudopotential description at $`p=0`$, the situation at high compressions will not be too different from that at $`p=0`$. As a result, one can expect that the information about the lattice properies of Ir (and its analogue Rh) at $`p0`$ which was derived from the pseudopotential calculations appears ot be raliable.

The investigation described in the present paper became possible in part thanks to the financial support of the International Science Foundation (grant RGQ300) and the Russian Foundation for Fundamental Research (grant 95-02-06426).

References

Browman, E. G., and Kagan, Yu. M., 1974, Uspekhi Fiz. Nauk, 112, 369.

Browman, E. G., Kagan, Yu. M., and Kholas, A., 1971, Zh. Eksp. Teor. Fiz., 61, 737.

Ceperly, D. M., and Alder, B. J., 1980, Phys. Rev. Lett., 45, 369.

Fomichev, S. V., Katsnelson, M. I., Koreshkov, V. G., and Trefilov, A. V., 1990, Phys. Stat. Sol. (b), 161, 153.

Geldart, D. J. W., and Taylor, R., 1970, Canad. J. Phys., 48, 155.

Gornostyrev, Yu. N., Katsnelson, M. I., Mikhin, A. G., Osetskii, Yu. N., and Trefilov, A. V., 1994, Fizika metall. Metallovede., 77(2), 79.

Greenberg, B. A., Katsnelson, M. I., Koreshkov, V. G., Osetskii, Yu. N., Peschanskikh, G. V., Trefilov, A. V., Shamanaev, Yu. F., and Yakovenkova L. I. , 1990, Phys. Stat. Sol. (b), 158, 441.

Ivanov, A. S., Katsnelson, M. I., Mikhin, A. G., Osetskii, Yu. N., Rumyantsev, A. Yu., Trefilov, A. V., Shamanaev, and Yu. F., Yakovenkova, L. I., 1994, Phil. Mag., B69, 1183.

Katsnelson, M. I., Naumov, I. I., Trefilov, A. V., Khlopkin, M. N., and Khromov, K. Yu., 1996, submitted to Phil. Mag. B .

Katsnelson, M. I., Peschanskikh, G. V., and Trefilov, A. V., 1990, Fiz. Tverd. Tela, 32, 570.

Korenovsky, N. L. and Polyakova, V. P., 1990, in Single crystals of metals Moscow, Nauka, in Russian, p72.

Perdew, J. R., and Zunger, A., 1981, Phys. Rev., B23, 5048.

Sigalas, M., and, Papaconstantopulos, D. A., 1994, Phys. Rev., B50, 7255.

Sigalas, M., Papaconstantopulos, D. A., and Bacalis, N. C., 1992, Phys. Rev., B45, 5777.

Vaks, V. G., Kravchuk, S. P., and Trefilov, A. V., 1977, Fiz. Tverd. Tela, 19, 1271.

Vaks, V. G., Zarochentsev, E. V., Kravchuk, S. P., Safronov, V. P., and Trefilov, A. V., 1978, Phys. Stat. Sol. (b), 85, 749.

List of tables to the paper by Katsnelson et al “Thermal expansion …”

Table 1. The values of the phonon frequencies $`\omega (𝐪)`$ and the microscopic Gruneisen parameters $`\gamma (𝐪)`$ in the simmetric points points of the Brillouin zone for Ir and Rh. $`\omega (𝐪)`$ in THz. A-H stands for the values calculated with the Animalu-Heine pseudopotential, I stands for the values calculated with the addition proposed in I.

Table 2. The contribution of the zero point vibration to energy and pressure for Ir and Rh. A-H stands for the values calculated with the Animalu-Heine pseudopotential, I stands for the values calculated with the addition proposed in I.

Figure captions to the paper by Katsnelson et al “Thermal expansion …”

Fig. 1 Microscopic Gruneisen parameters for Ir calculated for the pseudopotential from Ivanonv et al 1994 (dashed line) and from I (solid line).

Fig. 2 Microscopic Gruneisen parameters for Rh calculated for the pseudopotential from Ivanonv et al 1994 (dashed line) and from I (solid line).

Fig. 3 Temperature dependence of the linear coefficient of thermal expansion (a) and the macroscopic Gruneisen parameter (b) for Ir. Dashed line - calculation for the pseudopotential from Ivanov et al 1994; solid line — calculation for the pseudopotential from I. $``$ in fig. 3a — experiment from Korenovsky and Polyakova, 1990. $`T_{pl}=\mathrm{}\omega _{pl}=\mathrm{}(4\pi Ze^2/M\mathrm{\Omega }_0)^{1/2}`$, $`\omega _{pl}`$ is the ionic plasma frequency.

Fig. 4 Temperature dependence of the linear coefficient of thermal expansion (a) and the macroscopic Gruneisen parameter (b) for Rh. Dashed line — calculation for the pseudopotential from Ivanov et al 1994; solid line — calculation for the pseudopotential from I. $``$ in fig. 4a — experiment from Korenovsky and Polyakova, 1990. $`T_{pl}=\mathrm{}\omega _{pl}=\mathrm{}(4\pi Ze^2/M\mathrm{\Omega }_0)^{1/2}`$, $`\omega _{pl}`$ is the ionic plasma frequency.

Fig. 5 Equation of state for Ir at $`T=0`$. Solid line — calculation for the pseudopotential from I; $``$ab initio  calculations.

Fig. 6 Dependence of the phonon contribution to pressure upon volume at $`T=295K`$ and $`T=T_m`$ for Ir (solid line) and Rh (dashed line). The calculations were made for the pseudopotential from I.

Fig. 7 Temperature dependence of the phonon contribution to pressure for Ir (solid line) and Rh (dashed line). The calculations were made for the pseudopotential from I. $`T_{pl}=\mathrm{}\omega _{pl}=\mathrm{}(4\pi Ze^2/M\mathrm{\Omega }_0)^{1/2}`$, $`\omega _{pl}`$ is the ionic plasma frequency.

Fig. 8 Equation of state for Ir at T=295 K and $`T=T_m`$. $``$ is ab initio  calculation of $`p_0`$ with pseudopotential calculation of $`p_{ph}`$ using the addition proposed in I at temperature 295 K. $``$ is the same at the temperature $`T=T_m`$. Solid line — all quantities were calculated with the pseudopotential proposed in I, temperature is 295 K. Dashed line — the same at the temperature $`T=T_m`$.

Fig. 9 Dependence of the elastic moduli $`B_{ik}`$ on volume for Ir. The calculations were made for the pseudopotential from I.

Fig. 10 The density of electron states (a) and the electron energy spectrum (b) in Ir at equilibrium volume.

Fig. 11 Density of electron states in Ir at different pressures.

table 1

|  | $`\omega (𝐪)`$ | | | $`\gamma (𝐪)`$ | | |
| --- | --- | --- | --- | --- | --- | --- |
|  | (1,0,0) | (1,1/2,0) | (1/2,1/2,1/2) | (1,0,0) | (1,1/2,0) | (1/2,1/2,1/2) |
| Ir A-H | 8.06 | 7.00 | 8.11 | 2.39 | 2.27 | 2.33 |
|  | 5.65 | 5.85 | 4.03 | 2.18 | 1.83 | 1.60 |
| Ir I | 6.98 | 5.90 | 7.44 | 2.90 | 3.01 | 2.55 |
|  | 5.16 | 5.22 | 4.01 | 2.89 | 2.17 | 1.72 |
| Rh A-H | 10.03 | 8.63 | 9.97 | 2.17 | 2.10 | 2.15 |
|  | 7.04 | 7.10 | 4.91 | 1.99 | 1.68 | 1.44 |
| Rh I | 7.62 | 6.08 | 8.51 | 3.12 | 3.76 | 2.51 |
|  | 5.93 | 5.64 | 4.86 | 3.34 | 2.37 | 1.64 |

table 2

|  | Ir A-H | Ir I | Rh A-H | Rh I |
| --- | --- | --- | --- | --- |
| $`E_{zp}/E_0`$$`10^4`$ | 2.49 | 2.53 | 4.54 | 3.78 |
| $`P_{zp}/P_0`$$`10^3`$ | 2.63 | 3.1 | 3.77 | 5.25 |