Mario Saavedra-Torres

Frederik Tielens
email: frederik.tielens@upmc.fr

Juan C Santos
email: jsantos@unab.cl

Juan Carlos

Dibenzyl Disulfide Adsorption on Cu(111) Surface: A DFT Study

Keywords: Copper, dibenzyl disulfide (DBDS), DFT, adsorption

The adsorption of dibenzyl disulfide (DBDS) on a Cu(111) surface model was investigated by using density functional (DFT) calculations, considering energetic and electronic aspects. Several complexes were generated, where the bridge, hollow hcp, hollow fcc and top adsorption sites were considered. The results show that the Cu-S interaction guides the final complexes, and a secondary π-Cu weak interaction confers an extra stability. The complexes were grouped as physi-or chemi-sorption according to their adsorption energy applying a distortion decomposition model, with a preference by a double interaction of S with Cu (i.e. hollow hcp, and bridge sites). A degree of disulfide bond dissociation was observed in the complexes, being correlated with adsorption energies. From an electronic aspect, it was found that the electronic flow from copper to DBDS occurs in the most stables complexes, checked with charge analysis.

These results are agreed with experimental revelations of copper corrosion on power transformers.

Introduction

The adsorption of sulfur containing molecules on gold surfaces have attracted much attention and this especially since the discovery of the self-assembly of alkyl thiols on Au(111). [1][2][3] Thiols adsorption has also been studied on the other coinage metals [4][5][6][7] . Thiols 4,[8][9][10][11][12][13][14][15] , disulfides [16][17][18][19] , methionines [20][21][22][23][24] , cysteines [25][26][27][START_REF] Hoffling | Adsorption of Cysteine on the Au(110)-surface: A Density Functional Theory Study[END_REF][START_REF] Fischer | [END_REF][30][31][32][33] , etc. have been studied experimentally and theoretically. The reactivity towards sulfur containing molecules decreases going down the column of the coinage metals, explaining why molecular adsorption has been successfully observed in the perfectly ordered thiol chains on Au(111) (See reviews on the topic cited higher). On Cu(111), the more reactive surface of the three metals (Cu, Ag, Au), a different picture of the adsorption is expected and observed 5,34,35 . The surface itself will reconstruct and the adsorbed molecules easily dissociated.

Thus, understanding the mechanism of interaction of sulfur containing molecules is expected to be of high relevance for opening new perspectives towards improving the reactivity or stability of materials used in different applications 36 . Concerning the copper surfaces, sulfur containing molecules induce copper metal corrosion, which is a chemical phenomenon that triggers serious failures in industrial applications, and especially in power transformers [START_REF] Tumiatti | Conference Record of the 2008 IEEE International Symposium on Electrical Insulation[END_REF][START_REF] Lukic | [END_REF] .

Dibenzyl disulfide (DBDS) is an additive with antioxidant properties, frequently used in insulating mineral oil employed in power transformers. Under operation conditions, DBDS is known to further copper corrosion in electric equipment, by forming copper (I) sulfide Cu2S as the main product, with the subsequent production of other sulfur compounds such as benzyl sulfide and dibenzyl sulphide. In spite of successive studies about the chemical phenomenon, the mechanism remains unclear [39][START_REF] Ahmed Khan | International Conference on Advances in Power Conversion and Energy Technologies (APCET)[END_REF][START_REF] Oweimreen | [END_REF][42] . Since DBDS is related to the commonly used thiol molecules, in the well-known thiol self-assembled monolayers (SAM) on metal surfaces, 9,[43][44][45] , the analysis of the metal-sulfur bond is important to complete the already known information on the chemical properties of this type of systems.

In the present paper, the nature of the interaction between Cu(111) and DBDS is theoretically studied. A slab model is used to simulate the copper surface at reasonable accuracy, i.e. using periodic PBE-D2. In particular we clarify some aspects associated with copper corrosion by DBDS at the DFT level by accounting geometric, energetic and electronic properties that govern the interaction between DBDS and copper. The main focus of attention is to characterize locally the copper-sulfur interaction/bonding.

Computational Details

Computational Method

Calculations were performed in the frame of periodic DFT by means of the Vienna Ab Initio Simulation Package (VASP 5.3). 46,47 The electron-ion interactions were described by the projector augmented wave (PAW) method 48,49 , representing the valence electrons. The convergence of the plane-wave expansion was obtained with a cut-off of 400 eV. The Perdew-Burke-Ernzerhof (PBE) 50,51 , a Generalized Gradient Approximation (GGA) functional was used.

The adsorption of DBDS (See Fig. 1) on Cu(111) was modeled using a c7×7 periodic unit cell, with the aim to the study the adsorption of an isolated adsorbate molecule interacting with the surface. The Cu slab contains three layers, from which the bottom one was not allowed to relax and kept at the bulk positions. The lattice parameter was initially fixed at the experimental value, 3.61 Å, to build the initial slab 52 . Then, it was re-optimized and a subestimation of 3% from the original supercell was obtained. This small change is accepted because is in good agreement with the experimental parameter and the obtained theoretically from PBE 53 .

The sampling in the Brillouin zone was performed employing 4 and 13 k-points, resolved on 2×2×1 and 5×5×1 grids for the geometry optimizations (including dispersion corrections) and energy evaluation (at single PBE level), respectively.

Since our system involved organic molecules interacting through weak forces, the pure DFT energies obtained in periodic DFT PBE should be corrected. Therefore, we included Grimme D2 corrections 54 , which can be calculated using the presently used VASP version, although, D2 corrections overestimate the binding energy ( 55,56 and references therein).

The optimization procedure consisted of locate initially the DBDS molecule at 2.3 Å from the surface, at the beginning of the geometry relaxation. This distance was taken from chemi-and physisorption data reported in several studies involving the adsorption of sulfur molecules on copper surfaces 34,35,57,[START_REF] Crutcher | electric-connectionco[END_REF] .

Calculation of the adsorption energy and charge transfer

Each complex can be characterized by its adsorption energy, Eads, which is calculated from the total energy of the ground state optimized geometries of the complex, Cu(111) slab and DBDS as following:

∆Eads = ECu(111)-DBDS -[ ECu(111) + E(DBDS) ]

(1)

∆Eads.D2 = ∆Eads + [ ECu(111)-DBDS .D2 -[ ECu(111).D2 + E(DBDS).D2 ] ] (2)
Eq (2) includes the dispersion contribution D2, explicitly for each term.

However, this term can be scaled to obtain a better energetic estimation, avoiding the most part of the overestimation of empirical of the metal bulk. A common used correction is that the adsorption energy is approximated to include only the dispersion from the top layer of the surface and the adsorbate (DBDS), defining the term Eads.D2.1layer as:

∆Eads.D2.1layer = ∆Eads + [ (ECu(111)-DBDS.D2 ) -[ ECu(111).D2 + E(DBDS).D2] (3)

Moreover, an energetic decomposition scheme can be inspected, to describe energetically the change along the adsorption process, through a two part scheme of deformation/interaction [START_REF] Scaranto | [END_REF] .

This partition scheme consists of two contributions, deformation (∆Edef), and interaction energy (∆Eint), related as:

∆Eads = Edef + Eint (4)
Where ∆Edef ( X ) = E(Xdeformed )-E(Xequilibrium) ; X=Cu(111) , DBDS (5) Edef precise an energetic measure of perturbation over the gas phase equilibrium geometry of each molecule, to obtain their respective geometry in the complex. Is clear, that Eint, represents the interaction energy between the deformed molecules in the complex and can be easily computed, from eq. ( 4) and ( 5).

The amount of total electronic charge transferred between DBDS and the Cu surface, was quantified by a Global Charge Transfer descriptor (GCT), which corresponds to the sum of atomic charges ( ) on each molecule (A), i.e. the DBDS or the surface:

GCT = ∑ ; ∈ (6)
For this purpose the Bader Charge Analysis 60,61 was used.

Additionally, the Electron Localization Function Its numerical values are conveniently mapped on the interval (0,1) facilitating its analysis.

According to the interpretation of the ELF, a region of the space with a high value of ELF corresponds to a region where it is more probable to localize a pair of electrons of opposite spin.

The topological analysis of the ELF gradient field [63][64][65] provides a mathematical model enabling the partition of the molecular position space in a set of continuous and non-overlapping basins of attractors that present in principle a one-to-one correspondence with electron pairs.

In this way, an accurate calculation of chemical local objects such as bonds, lone pairs, or atomic shells can be achieved. The basins are either core basins labeled C(A) or valence basins V(A,...) belonging to the outermost shell. Valence basins are characterized by their coordination number (synaptic order) with core 66 . The original work of Silvi and Savin on the ELF generated a fruitful field of applications in a variety of chemical problems, ranging from structural and chemical reactivity studies as well as the study of chemical reactions [67][68][69][70] . This scheme will be applied to evaluate the S-S bond in the complex.

q A 6
Results and discussion

Geometric and Energetic Aspects

In this work, 20 optimized adsorption complexes were defined combining the different possible S adsorption sites on Cu(111) surface (top, hollow hcp, hollow fcc, and bridge) with the three isoenergetic DBDS gas phase conformations, obtained from a previous study 71 (See Fig. 1).

Using PBE, the structure B was more stable than Typically, the adsorption site for both S atoms is found close to a bridge site. Other adsorption sites slightly higher in energy are also present, such as on top (close). Different adsorption geometries (mainly due to different benzyl conformations) showing different adsorption energies can adsorb on the same adsorption site. In order to rationalize the data, we defined the vertical spacing distance z and the shortest S-Cu distance, d1 and d2 (See Fig. 2, Table 1).

The most stable physisorbed geometries, shown in Fig. 3, can be divided in two groups: a) having an S-Cu atop adsorption and b) 2S-Cu atop adsorption. The first is slightly more favorable (structure B7). It should be noted that the benzyl groups tend to orient parallel to the Cu(111) surface.

From Table 1, the S-S distances (d3) between 2.0 and 2.2 Å, corresponds to physisorbed complexes and also correlates with a high average values of d1 and d2 ( > 2.56 Å), with adsorption sites that vary between bridge, top and hollow sites. In the other side, the group of structures having d3 distance larger than 3.4 Å and an average of d1 and d2 lower than 2.56 describes the chemisorption complexes, where sulfur atom prefers the adsorption on bridge (slightly to hollow) sites.

In summary the most favorable adsorption sites was found to be bi-coordinated bridge in the case of physisorption and a bridge (slightly to hollow) for each dissociated fragment of benzylthiol, in the case of chemisorption (See Fig. 3 and Fig. 5).

The z distance is related to the strength of the interaction of the sulfur atom with the Cu(111)

surface. Moreover, the difference of the vertical spacing between the S atom and the closest Cu surface atom and z give a geometrical measure for the deformation of the surface. After adsorption, and especially with dissociative chemisorption, a Cu atom can be lifted up out of the surface, and can be considered as a precursor state to the formation of an adatom. Similar results have been reported in the adsorption of thiols on Au and Ag surfaces, being less evident than Cu, attributable to the metal hardness. Looking at z in Table 1, one can conclude that the effect of chemisorption on the z coordinate of the surface copper atoms is particularly noticeable over ∆z1, where values of ∆z1 lower than 5.98 corresponds to chemisorbed species, while higher than this value are related with physisorption complexes. The most affected copper atoms of the top layer, are lifted up by 0.2 Å. The magnitude of surface deformation due to the tension generated by the adsorption of DBDS will certainly increase with increasing coverage, as is noticed in thiol/Au(111) SAMs. 9,72 This surface tension is observed over a range of copper atoms in the neighborhood of the binding sulfur atom (See Fig. 4).

The adsorption energies (See Table 2) are in agreement with the geometrical parameters. Indeed the adsorption complexes having the shortest distances between the DBDS sulfur atom and the Cu(111) surface (d1,d2), correlates with the highest (absolute value) adsorption energies. The adsorption energy helps to identify that the phenyl-Cu proximity, shows a weak π-Cu interaction (see lower), leads to an extra stabilization of the adsorption complexes (structures B1, B2, and B3), which correspond to the chemisorbed species (See Fig. 5).

In general the complexes formed from conformer type B and C (See Fig. 1), have more stable adsorption energies due to the proximity of the phenyl groups to the surface, favoring -Cu interactions.

The deformation/interaction decomposition scheme displayed in Table 3, shows that in case of physisorbed complexes, the surface is slightly more perturbed than DBDS. The low values of ∆Edef and almost the total contribution of ∆Eint to ∆Eads, explain the low geometrical effect in the adsorption. In the other side for chemisorbed complexes, the ∆Edef are higher and particularly ∆Edef(DBDS) is greater than the surface and it represents more than 70% of ∆Eads with the inclusion of an important part of the dissociative process. Additionally, ∆Edef(DBDS) represents more than the 50% of absolute value of ∆Eint, while the ∆Edef (Cu) is less than 5% for ∆Eint.

Those assesses, explains an important role of the geometry relaxation of DBDS related with dissociation over the adsorption.

Electronic Aspects:

Global Charge Transfer

A global charge transfer descriptor was introduced to quantify the relation between the adsorption energies and charge transfer between DBDS and the Cu(111) surface. Table 2 summarizes this quantity. It can be seen that the direction of the charge transfer is highly related to the kind of adsorption (physi-or chemi-sorption). In the most stable cases (B1, B2, and B3, i.e. chemisorbed species), the DBDS molecule acts as an electron acceptor. Most transferred charge is then allocated on the disulfide group, reaching a maximum when dissociation of the S-S bond occurs.

Conversely, in the physisorption processes the DBDS molecule act as donor but with a lower rate of charge transfer to the Cu(111) surface (See Fig. 6).

Also, a spontaneous (barrier less) dissociation of the disulfide bond is observed for chemisorbed species, which increases the ability of DBDS to become an electron acceptor. In this case, the S-S distance increases from 2.0 to 2.2 Å and from 3.5 to 4.5 Å for physisorption and chemisorption, respectively. Finally, bi-coordinated interactions of the sulfur atoms with the Cu atoms and preferential orientation of the phenyl rings parallel to the Copper surface were observed.

Electron Localization Function (ELF) Analysis.

In the Fig 7 a representation of the ELF basins is shown to help visualize the most relevant regions of valence electron density that characterize the DBDS in its adsorption process on copper surface.

The basin population analysis using pseudopotentials is not accurate enough, mainly to describe the population associate to sulfur atoms. However, the topological description gives us important information about the bonds in the different obtained complexes in this work.

In all kind of formed complexes by physi-or chemi-sorption, each benzylsulfide moiety of the adsorbed DBDS presents two kinds of C-H bonds, aromatic (6) and methylene (2), represented by disynaptic basins V(H,C). The six carbon-carbon aromatic bonds are represented by disynaptic V(C-C), the shape and inclusive, the population close to 3e is characteristic of delocalized bond in aromatic systems, and the one non aromatic bond, V(C-C), has population close to 2.1 e.

Around the sulfur atom is possible to visualize the disynaptic V(C-S) with population between 1-1.5e and two monosynaptic basins V(S) which characterize the electron lone pairs. At least, one of this are located between sulfur and surface. But the main difference are centered in the disynaptic basin V(S-S), which describes the disulfide bond. In the physisorption complexes is possible to visualize it but in the chemisorption complexes this basin disappears. It indicates the disulfide dissociation in the most stable formed complexes.

Corroboration of topological description was made using the structure of isolated DBDS with the geometry adopted in the complex. First at the same level of calculation previously described and confirmed by all electron calculation using the same functional and def2-TZVP basis set. The same description was obtained, and the disynaptic V(S-S) appears and disappears in the DBDS physi-and chemi-sorption complexes, respectively.

Conclusions

The molecular (physisorption) and dissociative (chemisorption) adsorption of DBDS on Cu(111) was studied by means of periodic DFT. Physisorption was investigated in detail by optimizing systematically all possible conformations. The adsorption site was found to be bi-coordinated bridge for physisorption and on simple bridge for chemisorption. The adsorption energies were calculated and decomposed, validating the high contribution of relaxation in the adsorption process. The charge transfer and the ELF function were used to describe the electronic structure of the adsorption complex. It was found that a charge transfer from the Cu(111) surface to the disulfide group occurs, reaching a maximum in the structures where S-S bond dissociation is observed, whereas in the physisorption processes the DBDS molecule act as donor but with a lower rate of charge transfer to the Cu surface. The ELF shows the absence and presence of the S-S bond in chemisorption and physisorption, respectively.

The dissociative adsorption of DBDS on Cu(111) can be considered as the first step of the copper corrosion phenomenon in which an oxidation process takes place on the surface.

  (ELF) was studied to characterize the existence of the S-S bond while the DBDS interact with the surface. The ELF of Becke and Edgecombe 62 provides an orbital independent description of the electron localization based on strong physical arguments regarding the Fermi hole. The ELF is defined in terms of the excess of local kinetic energy density due to the Pauli Exclusion Principle and the Thomas-Fermi kinetic energy density.

  A and C conformers by 0.03 and 0.07 eV respectively. It confirms the isoenergetic character of three conformers. In the succeeding optimization procedures, the adsorption complexes between DBDS and Cu(111) results with different interaction sites, mainly noticed by S-Cu geometrical parameters. A schematic representation of showing selected geometrical parameters is presented in Fig 2.

  Fig 1. Three DBDS conformers of lowest energy used in this study as starting structure to generate the adsorption complex.

Fig 2 .Fig 3 .

 23 Fig 2. Schematic representation of the adsorption complex showing selected geometrical parameters used in the discussion of the adsorption site and geometry.

Fig 4 .Fig 5 .

 45 Fig 4. Adsorption geometry of the most stable dissociated DBDS molecule on Cu(111) showing the distortion through Cu-Cu distances (in red) after adsorption.

Fig 6 .

 6 Fig 6. General Charge transfer (GCT) vs adsorption energies (Eads and Eads_1layer_D2).

20 Fig 7 .

 207 Fig 7. Selected ELF Isosurfaces (ELF=0.80) for (A) Physi-and (B) Chemi-sorption complexes formed by DBDS adsorption on copper surface

Table 1 .

 1 Geometrical parameters describing the final adsorption geometries and sites: d1, d2, d3 and z (for sulfur 1 and sulfur 2, See Fig.2). (values in Å).

	Complex Final Adsorption Site (S1 S2)	d1	d2	d3 ∆z1 z2
	A1	bridge-top	bridge-top	2.75 2.98 2.09 6.47 6.29
	A2	bridge-top	bridge-top	2.78 2.66 2.09 6.19 6.28
	A3	top	top	2.96 2.99 2.09 6.48 6.39
	A4	bridge-top	bridge-top	2.97 3.02 2.09 6.51 6.42
	A5	hollow1-top	bridge-top	2.57 2.87 2.11 6.51 6.13
	A6	hollow1-top	bridge-top	2.56 2.95 2.14 6.48 6.07
	A7	bridge-top	bridge	2.87 3.06 2.08 6.71 6.42
	B1	hollow1-bridge hollow1-bridge 2.41 2.42 4.42 5.98 5.96
	B2	hollow1-bridge	hollow1	2.49 2.31 3.42 5.81 6.03
	B3	hollow1-top	hollow2	2.81 2.24 3.92 5.67 6.30
	B4	top	bridge-top	2.93 2.85 2.09 6.38 6.37
	B5	bridge-top	bridge-top	2.66 2.70 2.13 6.24 6.18
	B6	top	top	3.00 3.04 2.08 6.58 6.36
	B7	top	bridge-top	2.82 2.95 2.10 6.40 6.30
	B8	hollow2-top	hollow1-top	3.02 2.92 2.07 6.58 6.69
	C1	top	top	2.96 2.97 2.10 6.37 6.46
	C2	hollow2-bridge hollow1-bridge 2.58 2.61 2.18 6.14 6.09
	C3	top	bridge-top	2.95 2.76 2.09 6.29 6.43
	C4	hollow1-top	hollow2-top	2.77 2.77 2.15 6.40 6.42
	C5	hollow2-top	hollow1	3.00 3.05 2.07 6.76 6.65

Table 3 .

 3 Deformation energies (∆Edef) of surface and DBDS, and interaction energies (∆Eint) for complexes, in presence and absence of D2 vdW contribution (in eV). Chemisorption complexes are highlighted.

	Complex	∆Edef(Cu)	∆Edef(DBDS)	∆Edef (Cu+DBDS)	∆EintD2	∆Eint
	A1	0.15	0.07	0.22	-2.31	-0.66
	A2	0.15	0.06	0.21	-2.24	-0.50
	A3	0.11	0.11	0.22	-2.25	-0.67
	A4	0.15	0.02	0.17	-2.19	-0.69
	A5	0.17	0.11	0.28	-2.22	-0.57
	A6	0.18	0.09	0.27	-2.15	-0.56
	A7	0.14	0.04	0.18	-1.72	-0.36
	B1	0.21	3.64	3.85	-7.69	-5.02
	B2	0.28	3.13	3.41	-7.17	-4.39
	B3	0.35	3.45	3.80	-7.28	-4.68
	B4	0.15	0.20	0.35	-2.92	-0.48
	B5	0.18	0.27	0.45	-2.97	-0.29
	B6	0.17	0.26	0.43	-2.94	-0.24
	B7	0.14	0.18	0.33	-2.82	-0.40
	B8	0.08	0.10	0.18	-2.07	-0.03
	C1	0.15	0.31	0.46	-3.21	-0.32
	C2	0.17	0.31	0.48	-3.00	-0.33
	C3	0.14	0.16	0.30	-2.63	-0.52
	C4	0.20	0.09	0.29	-2.27	-0.34
	C5	0.08	0.04	0.12	-1.95	-0.18

Acknowledgements

The authors acknowledge the financial support by FONDECYT through the project number 1120785 and Universidad Andres Bello Grant DI-497-14/R. M.S.-T. thanks CONICYT for a Ph.D. scholarship and Universidad Andres Bello for support through grant DI40/12-I. F.T. thanks to FONDECYT for several appointments as Invited Professor (2013 and 2015).


Bernard Goossens

David Parello

Katarzyna Porada

Djallal Rahmoune

Parallel Locality and Parallelization Quality

Introduction

To run a program in parallel on a multi-core processor today, one must rewrite its C code to add by hand or through tools some OS parallelizing primitives such as pthread. Even with high level interfaces like OpenMP or MPI, parallelizing is not an easy job for two reasons: (i) if the resulting code is not enough synchronized, the computation is not deterministic and (ii) if it is too much synchronized, it is not parallel enough. This is illustrated by figure 1 [Copyright notice will appear here once 'preprint' option is removed.] showing the C implementation of a vector sum reduction and figure 2 showing its pthread implementation.

u n s i g n e d l o n g sum ( u n s i g n e d l o n g a r r a y [ ] ,

u n s i g n e d l o n g n ){ i f ( n==1) r e t u r n a r r a y [ 0 ] ; i f ( n==2) r e t u r n a r r a y [ 0 ] + a r r a y [START_REF] Baskaran | Automatic data movement and computation mapping for multi-level parallel architectures with explicitly managed memories[END_REF] ; r e t u r n sum ( a r r a y , n / 2 ) + sum(& a r r a y [ n / 2 ] , n-n / 2 ) ; } Figure 1: A vector sum reduction: C code The pthread coding is known to be tricky. A first (incorrect) version has no pthread join synchronization, which highlights (i). As a result, the run is not deterministic and the returned sum may be incorrect. A second version places the first pthread join synchronization on line 11 which serializes the second recursive call after the first one, which highlights (ii). As a result, the run is not parallel enough. Only the third version which places the two joins on lines 16 and 18 is satisfactory.

This example abstracts a first difficulty of hand-made parallelization with non deterministic OS primitives: achieve the exactly needed synchronization.

Edward Lee [START_REF] Lee | The problem with threads[END_REF] writes "Threads (...) make programs absurdly nondeterministic, and rely on programming style to constrain that nondeterminism to achieve deterministic aims.". Section 2 shows how to run deterministically in parallel.

Instead of hand-parallelizing the code, the developer may rely on its compiler to automatically do loop vectorization or loop parallelization. A simple example as the program given on figure 1 cannot be automatically parallelized by gcc. Irregular code structures are a second problem. Section 3 explains how loops are parallelized in our approach.

A third problem of parallelization is the memory organization of the data [START_REF] Baskaran | Automatic data movement and computation mapping for multi-level parallel architectures with explicitly managed memories[END_REF]. In the sum example, the array to be summed is declared for example as a global variable. Hence, it is centralized when the computation is distributed. As a result, each thread brings the array pieces it needs from DRAM, where it resides. Cache can help but neighbour cores are slowed down by memory contention and the cache miss rate is impacted by the array distribution. Moreover, in a program updating shared data, keeping caches coherent requires complex hardware which slows down average memory access time. Caches and memory hierarchy, as well as branch predictors, are hardware features that rely on the principle of locality, which in essence is founded on the centralisation of data (caches) and fetched code (predictors). When the code and the data are distributed, it is the parallel locality which applies to data. The parallel locality principle is that a consumer should be as close as possible from its producer. In this paper we propose a measure of the producer to consumer distance, which is a way to quantify the parallelization quality. In [START_REF] Goossens | Toward a core design to distribute an execution on a manycore processor[END_REF], a parallelizing core hardware is proposed to distribute an execution on a manycore processor. As the parallelization is dynamic, the three problems just mentioned are more easy to tackle. The programming model presented in this paper relies on this core hardware. The next section introduces its ISA 1 extension with the fork instruction, its new way of fetching in parallel to build a totally ordered trace and its parallelization of the renaming process, extended to memory locations.

A Deterministic and Parallel Run of C Code

Deterministic Parallel Execution

The sequential execution of the code on figure 1 is deterministic. The C code fixes a total order that the run follows. The core may run machine instructions out-of-order [START_REF] Tomasulo | An efficient algorithm for exploiting multiple arithmetic units[END_REF] [6], i.e. in a partial order derived from Read-After-Write (RAW) register dependences, but the total order semantic is preserved.

If a parallel execution is based on a total order, it is deterministic. A totally ordered trace can be built in parallel. For example in the C sum function, the trace can be built top-down. Such an outof-order construction of a totally ordered trace is possible because the control flow instructions are partially ordered. For example in the sum function, the control flow path in the second recursive call is independent from the control flow path in the first recursive call. This means that both calls traces can be built in parallel and orderly connected afterwards.

To avoid complications in the trace building, the hardware in [START_REF] Goossens | Toward a core design to distribute an execution on a manycore processor[END_REF] computes the control instructions targets rather than predicting them. Computing is slower than predicting but computing tens of branches in parallel is more efficient than predicting tens of 1 Instruction Set Architecture branches in sequence, parallelism being more cost-effective than a sequential predictor, even a perfect one.

The fork Machine Instruction

The first 18 lines on figure 3 are an x86 translation of the sum C code (AT&T syntax). The call and ret instructions are replaced by fork and endfork. The fork instruction semantic is to keep on fetching along the continuation path, i.e. go to the fork instruction label. Simultaneously, a new core starts fetching along the resume path, i.e. the instruction following the fork in the text. The endfork instruction semantic is to stop fetching along the current path. Contrarily to the call and ret pair, no return address is pushed/popped.

Moreover, the given code assumes the forked path (i.e. the resume code) receives a copy of the stack pointer (i.e. x86 register rsp), meaning that both paths use the same stack area. The hardware in [START_REF] Goossens | Toward a core design to distribute an execution on a manycore processor[END_REF] also copies rbp, rdi, rsi and rbx. These copies are better than push/pop because a push in a function prologue and a pop in its epilogue create RAW dependences between the epilogue and the prologue of the next function call, serializing them.

In the code on figure 3, register rsi and rdi hold arguments n and array. The sum function code leaves the computed sum in register rax (lines 3, 5 and 16), from where it is read by the resume path (lines 11 and 16).

Memory Renaming

In the hardware presented in [START_REF] Goossens | Toward a core design to distribute an execution on a manycore processor[END_REF] the trace is run in the partial order of its dependences. False register dependences are removed by renaming [START_REF] Sima | The design space of register renaming techniques[END_REF]. Many true register dependences are also removed by copying. For example, the computation of n -n/2 in line 13 (figure 3) depends on registers set in lines 7 and 8. As rbx and rsi are copied to the resume path, line 13 can be run in parallel with the continuation path in line 1.

False memory dependences are also removed. For example, as the stack pointer moves back and forth (allocation in line 10 and disallocation in line 17), different code portions use the same location, creating Write-After-Read and Write-After-Write memory dependences removed by renaming.

Memory renaming schemes have been proposed [START_REF] Moshovos | Dynamic speculation and synchronization of data dependences[END_REF][15], mainly to accelerate loads. The renaming relies on a predictor to quickly decide if a load depends on a previous store. However, a prediction based mechanism is not suited to eliminate false memory dependences. In [START_REF] Goossens | Toward a core design to distribute an execution on a manycore processor[END_REF], the memory hardware renaming is based on a search along the instruction trace total order.

Renaming is parallelized. The two recursive call destinations are simultaneously renamed. Destinations can be renamed in any order. Sources can be renamed out-of-order at the conditions that (i) the trace is totally ordered and (ii) all the destinations between a source and its producer are renamed.

The totally ordered trace is built from pieces which are fetched in parallel and later connected. When fetched, each instruction is renamed. A renamed place is allocated to hold the destination and for the sources, the closest producer is looked for by a backward search through the built trace. If a piece of the trace is missing, the search is suspended until the trace is extended. The fetch of next instructions continues during the source renaming search.

Parallel Construction of the Totally Ordered Trace

Figure 4 shows the parallel fetch of the sum function. The fetch is distributed on 11 cores 2 . It is done in 7 successive steps. At the first step, only core 1 is fetching. It fetches the start of the sum function code at line 1 (the instructions are the ones on figure 3). The rsi register holds the number n of elements to be summed, i.e. n = 10. Instructions on lines 1 and 2 are fetched and computed, requiring 2 They may be SMT-like logical cores. Instruction on line 9 (figure 3) is a fork. Core 6 3 (figure 4) receives the resume address from core 1 (i.e. line 10). At step 2, core 1 fetches the continuation path (line 1) and core 6 fetches the resume path (line 10). At step 3, four cores are fetching in parallel. Each core is able to compute its own control path and continue fetching independently from other cores.

= i + n * / movq %r b x , %r s i / * n = n -n /2 * / f o r k i n i t / * i n i t (& a r r a y [ n / 2 ] , n-n / 2 ) * / e n d f o r k / * s t o p * / main : movq $0 , %r b p / * i = 0 * / movq $10 , %r s i / * n = 10 * / movq $ a r r a y , %r d i / * r d i = a r r a y * / f o r k i n i t / * i n i t ( a r
v o i d f ( u n s i g n e d i n t i ){ a r r a y [ i ] = i ; } v o i d f o r r e c u r s i v e ( u n s i g n e d i n t f , u n s i g n e d i n t l , v o i d ( * body ) ( ) ) { i n t n = l -f +1; i f ( n == 1 ) { ( * body ) ( f ) ; r e t u r n ; } i f ( n == 2 ) { ( * body ) ( f ) ; ( * body ) ( f +1); r e t u r n ; } i f ( n != 0){ f o r
The path followed by a core is a section. A section starts after a fork instruction along the resume path and ends when an endfork instruction is reached. On figure 4, there are 11 sections, i.e. one per core. Sections may be long (e.g. a recursive descent like the section on core 1) or short (e.g. the transmission of the partial sum on core 4). On the average, sections are around 10 instructions long resulting in a fine grain parallelization.

The hardware in [START_REF] Goossens | Toward a core design to distribute an execution on a manycore processor[END_REF] builds the trace total order by linking the sections. Each section is linked to its successor and predecessor (dashed bidirectional lines on figure 4) 4 . Sections belonging to the same hierarchical level are also linked (plain lines, level predecessor). For example, cores 1, 6 and 11 fetch the highest level of the sum function, i.e. lines 1, 2, 7 to 9 (core 1), 10 to 15 (core 6) and 16 to 18 (core 11). The three involved sections are linked (plain lines). These direct links help finding stack renamings, bypassing the subtrees. For example, instruction 16 on core 11 finds the first half sum on top of the stack, saved by instruction 11 on core 6, without waiting for the construction of the cores 7 to 10 trace sub-tree.

A lot of instructions don't even need any source renaming as their sources have known values. In the sum example, core 1 does not rename any source as a copy of register rsi (i.e. rsi = 10) is received from the sum function caller. Only sources referencing a production of another section need to be renamed. For example, core 6 renames rax on line 11. This renaming is to be provided by core 5 which is core 6 closest predecessor producing rax on line 16. The predecessor link from core 6 to core 5 is established at step 4, when core 5 has reached its endfork instruction (line 18). Core 6 3 The core numbers are choosen arbitrarily 4 Successor links serve to export retired data sends a request to read rax to core 5 which returns the rax value to core 6 when it is computed.

In the sum example, there are 25 renamed sources. Ten of them concern array elements (lines 3 and 5 on figure 3). Five are references to the stack (i.e. (rsp), line 16). The last ten are references to register rax (lines 11 and 16).

The renamings of the array elements references deserve a particular treatment to avoid a long distance search of their producer in the trace. This is explained in section 4. The renamings of rax references need a round-trip communication between the renaming section and its predecessor (one way sends a request to read rax and the return way sends rax value). The renamings of the stack references 0(rsp) need a round-trip communication between the renaming section and its predecessor at the same hierarchical level (level predecessor). For example, line 16 in core 11 renames 0(rsp) from producing line 11 in core 6.

Full Renaming and Memory Management

As all destinations are renamed, the trace has a Dynamic Single Assignment form [START_REF] Vanbroekhoven | Transformation to dynamic single assignment using a simple data flow analysis[END_REF]. As described in [START_REF] Goossens | Toward a core design to distribute an execution on a manycore processor[END_REF], each core keeps its renamings in locally allocated storing resources. The cores have no data caches. The intermediate computations are kept in renaming storage until they are freed. Renaming storage is allocated up to the core capacity. When the core is full, its fetch is suspended 5 .

Only the final computations are saved to physical memory, when retired. As a result, the processor memory is naturally coherent as there is a single writer.

Such parallelizing cores are simpler than actual speculative cores as they need no branch predictor, no data cache and no vector computing unit. Hence, they are better suited to be the building block of future manycore processors.

Parallelizing Loops

Single for Loops

A single for loop is transformed into a divide-and-conquer tree of parallel iterations. For example, figure 5 is the transformation of the array initialization loop. Figure 6 shows the x86 code produced by the compiler for the for recursive function translating the for loop. Figure 7 shows the distributed run of the loop (only the 5 first iterations are fully shown; the 5 last ones -b5 to b9-are folded; the full tree has 9 steps distributed on 20 cores 6 ).

The body calls can communicate (i.e. body(j) can import a value from body(i) for all i < j). Each body section has its own control flow. A parallelizable loop has independent iterations. In this case, there is no communication between the body sections and they all have independent control flows.

Nested Loops

Nested for loops are transformed into two nested divide-andconquer trees of parallel iterations. For example, figures 8, 9 and 10 show the translation of two nested loops initializing a matrix. The x86 code is not shown but is easy to build. The run is fully parallelized with 100 sections organized as a binary tree. It is still easy for iterations of the inner loop to communicate. It is less easy for iterations of the outer loop as sections of intermediate inner loops may form a large separation of a consumer from its producer. In this case, the stack may help (the producer pushes and the consumer pops). Moreover, section 4 presents a general programming technique to optimize inter-sections communications.

While loops

Figure 11 shows a while loop. Figure 12 shows the transformation of the while loop into a for loop. The for loop is then transformed in a recursive function with a continuation condition (function cont cond). The run deploys a binary tree of l -f + 1 calls to the for break recursive function. The calls are run in parallel.

In [START_REF] Benabderrahmane | The polyhedral model is more widely applicable than you think[END_REF], the authors give many solutions to automatically transform loops that are difficult to parallelize at compile time.

A Programming Technique to Increase the

Quality of Parallelization

The Quality of Parallelization

A well parallelized execution should be composed of many parallel sections (fine grain is better than coarse grain), i.e. the average size sa of the sections should be low. This allows a simultaneous fetch on all the cores, i.e. fixing the Instruction Per Cycle (IPC) peak value to the number of cores. The communications between the sections should be rare and short distance. In other words, the average number of communications per section ca should be low and the average distance da from the producer to the consumer should be short, i.e. the number of visited sections should be low. Eventually, the number of instructions run ni should be low. Among these To increase the parallelization quality, we can decrease sa (i.e. either by increasing the number of sections or by decreasing ni). We can also decrease ca and da (increasing the number of sections should not increase the average communication distance). One way to keep da low is to recompute a value each time it is used. This should not increase ni too much (the recomputation should not be complex).

s t a t i c i n l i n e v o i d i n i t ( i n t i , i n t j ){ m a t r i x [ i ] [ j ]= i * 10+ j ; } s t a t i c i n l i n e v o i d empty ( i n t f , i n t l ){} main ( ) { u n s i g n e d l o n g i , j ; / * f o r ( i =0; i <10; i ++) f o r ( j =0; j <10; j ++) m a t r i x [ i ] [ j ]= i * 10+ j ; * / //
} } v o i d q u i c k s o r t ( i n t f , i n t l ){ . . . // f o r ( i 1 = f ; // i 1 < l && ( ( x1 = a r r a y [ i 1 ] ) <= p ) ; // i 1

An Example: the sum Function

Figure 13 shows a main function calling sum. Figure 3 is its x86 translation. The for loop is translated into a parallelized function init. The init function fetched trace has the same sections decomposition as the sum function. Figure 4 can be adapted to the init function fetch. For example, the upper leftmost box corresponds to the fetch of lines 19, 20 and then from 26 to 28.

Figure 14 shows the section decomposition of the main function run. A starting section (upper left rectangle) is continued by the figure 4 upper leftmost rectangle (hence, lines 35 to 38 are followed by the init function tree on figure 4 and all belong to the first section). This section forks a new section to start the execution u n s i g n e d l o n g n , l o n g ( * g e t ) ( ) , u n s i g n e d l o n g i ){ i f ( n==1) r e t u r n ( * g e t ) ( a r r a y , i ) ; i f ( n==2) r e t u r n ( * g e t ) ( a r r a y , i ) + ( * g e t ) ( a r r a y , i +1); r e t u r n sum ( a r r a y , n / 2 , g e t , i ) + sum(& a r r a y [ n / 2 ] , n-n / 2 , g e t , i+n / 2 ) ; } s t a t i c i n l i n e l o n g s e t a r r a y e l e m ( l o n g a r r a y [ ] ,

u n s i g n e d l o n g i ){ a r r a y [ i ]= i ; } s t a t i c i n l i n e l o n g g e t a r r a y e l e m ( l o n g a r r a y [ ] ,

u n s i g n e d l o n g i ){ s e t a r r a y e l e m ( a r r a y , i ) ; r e t u r n i ; } v o i d main ( ) { l o n g s=sum ( a r r a y , SIZE , g e t a r r a y e l e m , 0 ) ; p r i n t f ( " s=%l d \n" , s ) ; } Figure 15: Modified main and sum functions with short producer to consumer distance of the sum function (hence, the sum is computed in parallel with the array initialization). The forked section starts at line 39 and continues as the sum function tree on figure 4. Hence, lines 39, 1, 2, from 7 to 9, 1, 2, from 7 to 9 and from 1 to 6 belong to the same section computing array[0]+array [START_REF] Baskaran | Automatic data movement and computation mapping for multi-level parallel architectures with explicitly managed memories[END_REF].

There are 105 instructions run for the sum computation and 93 for the array initialization, hence ni = 198 (neglecting the printf run). There are 23 sections, i.e. sa = 198/23 = 8.61 (less than 9 instructions per section on the average).

There are 90 register and memory sources in the init run 7 , all of them having a distance of 1 (i.e. requiring no communication).

There are 122 sources in the sum run, 15 of them having a distance of 2, 10 having a distance 8 of 12 and the remaining 97 having a distance of 1. Hence the average distance is da = 337/212 = 1.59. The average number of communications per section is ca = 25/23 = 1.09 (only 25 sources among 337 need an import from another section, i.e. have a distance greater than 1).

Modifying the sum Function to Enhance the Quality of the Parallelization

To get the array values, function sum sections send renaming requests which travel along the trace, visiting one core per section in the trip, i.e. requiring 12 core-to-core communications. If these distances can be shortened, i.e. if the sum function section consumer can be closer to its init function section producer, the execution time will be reduced. Reducing da enhances the parallelization quality.

Figure 15 shows modified main and sum functions. A value is recomputed each time it is referenced rather than being transmitted. In manycore processors, local computation is cheaper than distant communication.

The array initialization is fused in the sum function. The sum function gets each array element. The get function calls a set function which sets the array element. Hence, each time an array element is read, it is also written. This ensures that a consumer and a producer belong to the same section.

As a general rule, a function has new arguments: a pointer on a get function to read variables used in the computation, a structure to encapsulate the arguments of the get function, a pointer on a put function to use computed variables and another structure to encapsulate the arguments of the put function.

Figure 16 shows the x86 code translation of the program depicted on figure 15.

Figure 17 shows the run trace of the modified sum program. The number of instructions run is ni = 122 (38% less work). There are 11 sections, i.e. sa = 122/11 = 11.09. The number of sections is reduced and the average size of sections is increased (i.e. the run is artificially less parallel because the reduction of the work has reduced the number of sections). There are 146 sources, among which 15 have a distance 2 and the others have all a distance 1. The average communication distance is da = 161/146 = 1.1. The number of communications per section is ca = 15/11 = 1.36.

The key factor is da, 31% reduced. All sections get their sources locally or from the preceding section, requiring a round-trip coreto-core communication.

Instead of having two separate computations, one to initialize the array and one to compute the sum of its elements, we have a single one, fusing the initialization into the computation. The programming technique which favours parallelization mimics a dataflow style: variables are not initialized and later used but their initialization is delayed until they are used.

When a variable is to be used multiple times, it can either be stored or recomputed. It is the compiler's job to estimate if the recomputation is better (i.e. reduces da enough without increasing ni too much). To be accessed fastly, a scalar variable can be stored in the stack. To look for it, a renaming request travels along the level predecessor links, bypassing sub-trees. 7 Including implicit sources like register ef lags. 8 Line 3 reads array[i] in the sum function from the value written by line 21 in the init function. The reading section is in the sum tree and the writing section is in the init tree. There are 12 sections from the writing one to the reading one.

Improving a Matrix Multiplication Program

Figure 18 presents a non optimized C code to compute matrix multiplication. Optimized versions take advantage of cache locality to reduce the average access time to the input matrix elements. In a parallel environment, very distant elements are simultaneously requested (or with a short time gap) and in this case, the cache does not help. Instead, it is a better policy to keep each producer close to its consumer.

When Figures 19 and 20 are the modified implementation to reduce da, the distance from producer to consumer.

The first part of the code defines two types to encapsulate the needed values for the get and put functions. It also contains Instead of an initialization in the matrix declarations, the code provides special functions. From these, it is possible to produce in a few instructions any input matrix element in any place of the code (the set mat a elem function fixes an element of matrix a in 4 instructions, using a table of branches issued from the switch control structures). The get mat elem function includes a call to the set mat a elem or set mat b elem function according to the referenced matrix.

The first part also contains the sum function, used to compute the sum part of the vector product of a line from matrix a with a column of matrix b. It is the same divide-and-conquer sum function as the one presented above, adapted to get the elements to be summed. The transmitted get function is get vec elem which reads a[i][k] and b[k][j] and return their products. Each section computing a product gets one element of matrix a and one element of matrix b from the same section running the set mat a elem or set mat b elem function. Hence, the distance da is no more a function of m, n and p but constant 1 for all these elements reads.

The modified example has an average distance which is ten times better than the basic one. However, each read requires 4 instructions instead of a single load in the basic trace. Hence, the number of instructions run ni is roughly 4 times bigger.

The second part of the program is the matrix multiplication function which sets each element of the result matrix c with the vector product computed by the sum function. The element set is transmitted to the put function, i.e. the print mat elem function which prints it in a formatted manner (n elements per line).

Related Works, Discussion and Conclusion

Automatic parallelization parallelizes as regular as possible loops with techniques based on the polyhedral model [START_REF] Feautrier | Dataflow analysis of array and scalar references[END_REF].

Parallelization based on OS threads suffer from the rather high overhead of OS primitives. Using gdb, we have measured the cost

Figure 2 :

 2 Figure 2: A vector sum reduction: pthread code

Figure 3 :Figure 4 :

 34 Figure 3: A vector sum reduction (x86 code)

  Figure 5: A parallelized for loop

  f o r r e c u r s i v e : movq %r s i , %r c x / * r c x = l * / su bq %r d i , %r c x / * r c x = l -f * / addq $1 , %r c x / * r c x = l -f +1 *

Figure 7 :

 7 Figure 7: Trace of for recursive(0, 9, body)

Figure 8 :Figure 9 :

 89 Figure 8: Inner for loop

Figure 10 :Figure 11 :

 1011 Figure 10: Two nested for loops

Figure 12 :

 12 Figure 12: The while loop is transformed into a for loop with a continuation condition

Figure 13 :Figure 14 :

 1314 Figure 13: A main function calling the sum function

Figure 16 :

 16 Figure 16: The x86 translation of the modified sum program

Figure 17 :

 17 Figure 17: Trace of sum(array, 10) (modified sum program)

  Figure18presents a non optimized C code to compute matrix multiplication. Optimized versions take advantage of cache locality to reduce the average access time to the input matrix elements. In a parallel environment, very distant elements are simultaneously requested (or with a short time gap) and in this case, the cache does not help. Instead, it is a better policy to keep each producer close to its consumer.Whenc[m][n] = a[m][p] * b[p][n] is parallelized, the run creates m * n * p sections, each computing a a[i][k] * b[k][j] product.Each product computation reads two elements which are looked for in the trace. Each element is set in the OS process start code which copies the values of the input matrices from the executable file to memory (it is the linker translation of lines 23 and 24 to initialize matrices

  Figure 18: A C program to multiply matrices

  #i n c l u d e <s t d i o . h> 2 v o i d i m a t m u l ( i n t * a , i n t * b , i n t * c , 3 u n s i g n e d i n t m, u n s i g n e d i n t n , u n s i g n e d i n t p ){ 4 // c [m ] [ n ] = a [m ] [ p ] * b [ p ] [ n ] 5 u n s i g n e d i n t i , j , k ;

	6	f o r ( i =0; i <m; i ++)
	7	f o r ( j =0; j <n ; j ++){
	8	* ( c+i * n+j ) = 0 ;
	9	f o r ( k =0; k<p ; k++)
	10	* ( c+i * n+j ) += * ( a+i * p+k ) * * ( b+k * n+j ) ;
	11	}
	12 }
	13	v o i d
	14	p r i n t m a t ( i n t * a , u n s i g n e d i n t m, u n s i g n e d i n t n ){
	15	u n s i g n e d i n t i , j ;
	16	f o r ( i =0; i <m; i ++){
	17	f o r ( j =0; j <n ; j ++)
	18	p r i n t f ( "%d " , * ( a+i * n+j ) ) ;
	19	p r i n t f ( "\n" ) ;
	20	}
	21 }
	22	main ( ) {
	23	i n t a [ 2 ]

8 2015/11/19 1

2015/11/19

As in actual out-of-order cores when no more renaming register is free.

Among which 10 do nothing more than fetching an endfork instruction. The compiler can eliminate such empty sections by reversing the continuation and resume paths. 4 2015/11/19

2015/11/19

} } s t a t i c i n l i n e i n t g e t m a t e l e m ( S MAT ELEM * s ){ i f ( s->a==&a [ 0 ] [ 0 ] ) s e t m a t a e l e m ( s->i , s->j ) ; i f ( s->a==&b [ 0 ] [ 0 ] ) s e t m a t b e l e m ( s->i , s->j ) ; r e t u r n * ( s->a + s->i * s->n + s->j ) ; } s t a t i c i n l i n e i n t g e t v e c e l e m ( S VEC ELEM * s ){ S MAT ELEM s t ; i n t l , r ; s t . a=s->a ; s t . i=s->i ; s t . j=s->k ; s t . n=s->p ; l=g e t m a t e l e m (& s t ) ; s t . a=s->b ; s t . i=s->k ; s t . j=s->j ; s t . n=s->n ; r=g e t m a t e l e m (& s t ) ; * ( s->v)= l * r ; r e t u r n * ( s->v ) ; } i n t sum ( i n t a r r a y [ ] , u n s i g n e d i n t n , i n t ( * g e t ) ( ) , S VEC ELEM * s t , u n s i g n e d i n t i ){ i f ( n==1){ s t ->v=&( a r r a y [ 0 ] ) ; s t ->k=i ; r e t u r n ( * g e t ) ( s t ) ; } i f ( n==2){ i n t l , r ; s t ->v=&( a r r a y [ 0 ] ) ; s t ->k=i ; l =( * g e t ) ( s t ) ; s t ->v=&( a r r a y [START_REF] Baskaran | Automatic data movement and computation mapping for multi-level parallel architectures with explicitly managed memories[END_REF] ) ; s t ->k=i +1; r =( * g e t ) ( s t ) ; r e t u r n l+r ; } r e t u r n sum ( a r r a y , n / 2 , g e t , s t , 0 ) + sum(& a r r a y [ n / 2 ] , n-n / 2 , g e t , s t , n / 2 ) ; } v o i d i m a t m u l ( i n t * a , i n t * b , i n t * c , u n s i g n e d i n t m, u n s i g n e d i n t n , u n s i g n e d i n t p , On an x86 processor (ubuntu 14.04), we obtained respectively 726 and 565 instructions, leading for the sum program example to an overhead of 1291 * log n instructions for an array of n elements. The threaded version needs to sum at least 2000 terms to be costeffective against the sequential code. On the other hand, today's hardware make it possible to send a few bytes core-to-core in a few cycles (e.g. 3 cycles -send + route + receive-for two neighbour cores on a grid with a NoC).

The existing contributions [START_REF] Kim | An instruction set and microarchitecture for instruction level distributed processing[END_REF][9] [START_REF] Ranjan | Fg-stp: Finegrain single thread partitioning on multicores[END_REF][12] on a hardware approach to automatize parallelization suffer from the low basic Instruction Level Parallelism (ILP). The hardware based parallelization in [START_REF] Goossens | Toward a core design to distribute an execution on a manycore processor[END_REF] overcomes this limitation in 3 ways: (i) very distant ILP is caught because fetch is parallelized, (ii) all false dependences are removed through full renaming and (iii) many true dependences are removed by copying values. The remaining dependences in a run are true ones related to the sequentialities of the algorithm which the program implements. In such conditions, the authors in [START_REF] Goossens | Limits of instruction-level parallelism capture[END_REF] have reached a high ILP (thousands), increasing with the data size, on benchmarks issued from the PBBS suite (parallel applications; available at URL http://www.cs.cmu.edu/ pbbs/index.html).

The number of transistors on a chip allows the integration of thousands of simple cores. It is urgently needed that any program, including the OS itself, be parallelized. Parallelization should be done fastly and reliably, with reproducible computations, which is ensured if the run is deterministic.


Gustavo Mata Osoro

Damien Bregiroux
email: damien.bregiroux@upmc.fr

Mai Pham Thi

Franck Levassort

Structural and piezoelectric properties evolution induced by cobalt doping and cobalt/niobium codoping in BaTiO 3

Keywords: Perovskite, barium titanate, cobalt doping, piezoelectricity

Lead-zirconate-titanate PbZr (1-x) Ti x O 3 (PZT), materials for piezoelectric applications are the dominant piezoceramics since their properties are high. However PZTs need to be substituted by safer materials due to health and environmental problems of lead. In this work, the crystal structure modifications of BaTiO 3 induced by cobalt doping and niobium/cobalt co-doping were investigated in the range 0.1 -25 at.%. Low Co and Co/Nb doping level significantly improve piezoelectric properties of BaTiO 3 . The piezoelectric constant, obtained for BaTiO 3 doped with 0.075 at.% of Cobalt (d 33 = 200 pC/N), is twice of that measured for pure BaTiO 3 .

At high doping concentration, the structure of BaTi 1-x Co x O 3 and BaTi 1-x Co x Nb x O 3 change from tetragonal to hexagonal while BaTi 1-x Co x Nb 2x O 3 remains tetragonal and gradually turns into cubic. It suggests that the ionic charge defect is balanced by oxygen vacancies formation when Ti 4+ are substituted by Co 2+ /Co 3+ .

Introduction

Barium titanate BaTiO 3 (BT) was one of the first materials used for the fabrication of piezoelectric devices. Nevertheless, it was quickly replaced following the discovery of lead zirconate titanate (PZT) based piezoelectric ceramics. Since several years, BT ceramics have attracted much attention again with the aim of developing lead-free piezoelectric ceramics [1]. Since BT ceramics show intrinsically lower piezoelectric performances than PZT, it is of prime importance to enhance properties via the use of dopants or innovative processing [2,3]. Many dopants have already been tested, including isovalent [4], donor or acceptor dopants [5,6]. Doping effect on structural and ferroelectric properties of BaTiO 3 results to be very complex and is not yet well understood since most of time, several parameters change at the same time (nature and amount of point defects, grain size, crystal distortions…) [7]. Doping barium titanate with cobalt has been proved to lower the temperature of the ferro-paraelectric phase transition (Tc). Cobalt was also found to be an effective dopant to flatten the dielectric constant at low temperature [8]. Cui et al. showed that the relative dielectric constant at Curie temperature of Codoped BaTiO 3 is higher than that of the undoped sample [9]. Cobalt was used for its magnetic properties to prepare BaTiO 3 -based multiferroic materials [10]. In addition, cobalt seems also to have an effect on the processing of BT ceramics since the grain boundary mobility decreases with Co doping, due to the segregation at the grain boundaries of defects induced by doping [11]. On the other hand, pentavalent ion doped, such as Nb 5+ on the Ti 4+ site, could changes the normally insulating behavior of BaTiO 3 to semiconducting [12].

In this work, the crystal structure modifications of BaTiO 3 induced by cobalt doping and cobalt/niobium co-doping were investigated. Associated ferroelectric and piezoelectric properties were also characterized for some compositions.

Experimental

Compounds were synthesized at 1100 °C for 2 h by conventional solid-state reaction of BaCO 3 , TiO 2 , Co 3 O 4 and Nb 2 O 5 . The reacted powder was milled and uniaxially pressed at 60 MPa before being sintered in a platinum crucible at 1350 °C for 4 h in air. Bulk density was calculated from the Archimedes method in water. X-ray diffraction measurements were taken in a Panalytical X'Pert Pro diffractometer (45 kV, 40 mA, λK α1 (Cu) = 1.5405 Å/forward Ge (111)). Cell parameters of samples were determined by Rietveld refinement in Lebail mode using Fullprof software [START_REF] Rodriguez-Carvajal | 2k: Rietveld, Profile Matching and Integrated Intensity Refinement of X-ray and Neutron Data, V 1.9c[END_REF]. Samples were poled for electrical measurements at 150 °C by field cooling process at 600 V.mm -1 . Temperature dependence of the dielectric constant was measured at 1 kHz using a HP 4192A impedance analyzer, between 20 and 177 °C. Piezoelectric constant d 33 of the samples were measured by means of quasistatic d 33 meter (ZJ-3A) based on Berlincourt method. The electromechanical properties of the pellets were deduced from the measurement of the complex electrical impedance around the fundamental thicknessmode resonance. A HP4294A spectrum analyzer (Agilent Technologies Inc.) and an impedance test kit were used as the experimental setup. The theoretical behavior of the electrical impedance was computed as a function of frequency for the radial and thickness-modes using a one dimensional model based on the Krimholtz-Leedom-Matthaei (KLM) equivalent electrical circuit [START_REF] Krimholtz | [END_REF]15]. A fitting process was used to deduce successively the radial and thickness-mode parameters of the pellets from the experimental data (sample's thickness << diameter) [16]. The radial-mode parameter k p (planar coupling factor) was obtained and the following thickness-mode parameters were deduced: the dielectric constant at constant strain ε 33 S /ε 0 , the thickness coupling factor k t , and the electrical loss factors tan δ e .

Results and discussion

Structural changes induced by Co doping and Co/Nb co-doping

Considering the valencies and ionic radii of Co, Ti and Ba, one can make the assumption that Co can be introduced in the ABO 3 perovskite structure on the B site. Therefore, the cobalt-doped BaTiO 3 powders were synthesized according to the following reaction:

BaCO 3 + (1-x)TiO 2 + x/3Co 3 O 4 → BaTi 1-x Co x O 3-δ + CO 2 + yO 2 (1)
The amounts of oxygen vacancies δ and dioxygen release y are dependant of the cobalt valence state in the final product. The formation of oxygen vacancies can be described according to the following reactions: tetragonal phase, h: hexagonal phase).

For Co +II : 2BaO + CoO → 2Ba Ba + Co Ti '' + 3O O + V O ∘∘ (2) For Co +III : 2BaO + Co 2 O 3 → 2Ba Ba + 2Co Ti ' + 5O O + 1 2 V O ∘∘ (3)
Fig. 1 shows a zoom on XRD patterns and the associated cell parameters of the synthesized powders for different values of x. Only two phases were detected, a tetragonal (P4mm, Z = 1) and hexagonal (P6 3 /mmc, Z = 6) perosvkite named t and h, respectively. The amount of the hexagonal form increases with cobalt doping, while that of the tetragonal one decreases concomitantly. For both phases, cell parameters remain within the error bar. Observed for the first time in 1883 and described in details in 1948 [17], the hexagonal phase of BaTiO 3 is stable above 1430 °C for undoped samples [18], but can be stabilized at room temperature by reduction of a part of Ti +IV to Ti +III [19] or by substituting Ti by dblocks elements (Mn, Fe, Co…) [20,21]. In contrast of the cubic or tetragonal form of the pervoskite structure in which all the TiO 6 octoedra are linked by corners, in the hexagonal form, 2/3 of the TiO 6 octahedra share a face to form Ti 2 O 9 dimers (Fig. 2). In these compounds, the presence of oxygen x in BaTi The chemical compositions of both phases are rather difficult to establish, but the negligible cell parameters variations of the tetragonal phase and the fact that the hexagonal one can be observed at very low doping level indicate that a major portion of the cobalt incorporates the hexagonal phase.

The existence of the hexagonal form of BaTiO 3 seems to be due to the presence of oxygen vacancies.

This phase should therefore not be observed if the charge compensation is made by codoping with a donor element like niobium on the Ti site.

The amount of Nb +V required for the charge compensation is dependent of the cobalt valence state, according to the following reactions:

For Co +III , Nb/Co = 1:

4BaO + Co 2 O 3 + Nb 2 O 5 → 4Ba Ba + 2Co Ti ' + 2Nb Ti ∘ + 12O O (4)
For Co +II , Nb/Co = 2:

b" a" c" 5BaO + 2CoO + 2Nb 2 O 5 → 5Ba Ba + 2Co Ti '' + 4Nb Ti ∘ + 17O O (5)
Fig. 3 shows the structural changes induced by Co/Nb codoping of BaTiO 3 , observed by XRD, for the two considered valence states of cobalt. When Nb/Co = 1 (Co +III ), the hexagonal form of BaTiO 3 can be observed since x = 0.005, suggesting the formation of oxygen vacancies in BaTiO 3 , probably because there is not enough Nb to fully balance the charge defects induced by the substitution of Ti by Co. In contrast, when Nb/Co = 2 (Co +II ), hexagonal form was never observed. Nevertheless, the structure of BaTiO 3 changes gradually from the tetragonal form to the cubic one. The broad peaks observed for intermediate values of x could be explained by composition heterogeneities within the samples.

Conclusions

The crystal structure modifications of BaTiO 3 induced by cobalt doping and niobium/cobalt co-doping were investigated. In Co doped BT samples, oxygen vacancies are responsible of the formation of the hexagonal form of BT. BT remains tetragonal or cubic for samples doped with 2 Nb for 1 Co. These results suggest that cobalt valence state in BT is +II. Low Co and Co/Nb doping level significantly improve piezoelectric properties of BaTiO 3 , and lead to a slight decrease of Tc.

Fig. 1 .

 1 Fig. 1. XRD patterns and associated cell parameters of BaTi 1-x Co x O 3-δ compounds as a function of x (t:

  key factor of the existence of the hexagonal form of the perovskite structure. It is well established that structure consisting of polyhedra with corner connection like tetragonal perovskite cannot accept high amount of oxygen vacancies. This can be explained by the strong coulombic repulsion between Ti +IV cations that should be induced by the absence of the bridging oxygen.Keith et al. showed that in the hexagonal form oxygen vacancies are located between the face sharing TiO 6 octahedra[22].

Fig. 2 .

 2 Fig. 2. Crystal structure of tetragonal (t, left) and hexagonal (h, right) perovskite.

Fig. 3 .

 3 Fig. 3. XRD patterns of Co/Nb substituted barium titanate.

Table 1

 1 These results suggest that, under the process conditions used in this work, the valence state of cobalt in BaTiO 3 is +II. These results are in accordance with those obtained by Aguilar who characterized by EPR technique the valence state of cobalt in Co doped BaTiO 3 single crystals and by Yu et al.[21,23] but disagree those recently obtained by Das and Roul[10].Main piezoelectric properties of three samples, i.e. pure, cobalt doped and cobalt/niobium co-doped BaTiO 3 , are reported in Table1. Because of very low resistivity and too high electrical losses, piezoelectric properties were not possible to be characterized for sample with cobalt doping level over than 0.01. Results clearly show that piezoelectric properties are significantly enhanced when BaTiO 3 is doped by cobalt. On the other hand, permittivity is lower and dielectric losses increase. The latter is consistent with the presence of oxygen vacancies[24]. Compared to cobalt doping, codoping with niobium does not seem to affect piezoelectric properties, but allows maintaining a good level of permittivity with low dielectric losses. As previously observed by Weill et al.[25], Curie temperature slightly decreases when BaTiO 3 is doped or codoped.Ferroelectric and piezoelectric properties of pure, Co doped and Co/Nb codoped BaTiO 3 ceramics

	3.2. Properties

τ: relative density, G: mean grain size, k p : planar coupling factor, ε 33 S /ε 0 : dielectric constant at constant strain; tan δ e : dielectric losses , d 33 : piezoelectric coefficient, k t : thickness coupling factor; Tc: Curie Temperature

Acknowledgments

This work was supported by the French Research Agency (ANR Hypercampus Mat&Pro 2010)