H Aknouche

A Outahyon

C Nouveau

R Marchal

A Zerizer

J C Butaud

JMATPROTEC Tool wear effect on cutting forces: In routing process of Aleppo pine wood

Keywords:

This paper uses the cutting forces in a routing process of Aleppo pine wood to estimate the tool wear effect. The aim is to obtain further information about the tool wear effect by monitoring the variation in the cutting forces. A Kistler 9257A 3 axes Dynamometer was positioned under the workpiece to measure the cutting forces at frequencies up to 10,000 Hz. The experiments were carried out on a CNC routing machine RECORD1 of SCM.

A carbide tool was used and the cutting parameters were fixed. The cutting speed was approximately 25 m/s. Dasylab software was used to capture the data. The results show a correlation between the tool wear and the computed angle (), between the tangential and cutting forces. In fact, the variation of () is unstable in the running period and stable in the linear wear zone, included in the interval [-1.11 • ; -1.10 • ]. This study was performed as part of a development program for the Algerian wood industry, hence the selection Aleppo pine wood as the working material.

Introduction

Since the first works of F.W. Taylor (1850Taylor ( -1915) ) a great deal has been learnt about tool wear, as it is a widely studied phenomenon. However, the development of new kinds of tools and new materials has expanded the experimental field. This has resulted in a renewed interest in the study of the phenomena of wear. This renewed interest is highlighted by examples of tool wear investigations in different application fields: [START_REF] Weinert | Considerations of tool wear mechanism when machining metal matrix composites (MMC)[END_REF] found, using new advanced materials like metal matrix composites, important differences in tool wear if the hardness of the reinforcement material is higher or lower than the hardness of the cutting material. In the case of dry machining of aerospace aluminium alloys, [START_REF] Nouari | Experimental analysis and optimisation of tool wear in dry machining of aluminium alloys[END_REF] show that the temperature is the major parameter inducing tool wear.

According to [START_REF] Roumesy | Remarks on the comparative representivity of wear law models of carbide tools[END_REF], the study of wear must be tackled by a gradual and systematic thought process i.e. first, systematically identifying the multiple manifestations of wear and observe their respective development. Then, determine the acceptable limits, according to the quality requirements. The limits will be in the form of quantitative criteria (e.g. crater wear and flank wear), which allow various manifestations of the macroscopic wear.

However, [START_REF] Anselmetti | Surveillance de la coupe et CFAO[END_REF] shows, depending on the various fields of application, that because of the relative nature of wear criteria, it is important to specify which are the most restrictive criteria. It is clear that, in a given field, the wear criteria chosen alone should be used to judge reliably when the tool must be removed from service for renewal of the cutting edge . The wear can be defined as the loss of material from the cutting edge due to mechanical or chemical factors associated with the cutting process. The cutting edge wear process according to the cutting path in machined material shows three areas that characterize the behaviour of the tool edge recession, which are; running (abrupt wear), linear wear (known as stability period) and catastrophic wear increasing significantly the recession of the tool edge (tool failure). [START_REF] Itaya | Development of CrN coated tool to create smooth surface of solid wood aiming to mit sanding process in furniture production[END_REF] concluded that the direct consequence of the cutting edge wear is the gradual loss of its ability to cut the machined material. The result being an increase in the cutting forces and the machine power consumption. These conditions also result in a poor wood surface quality. The wear of the cutting edge therefore has a big impact on the cost of the wood machining process.

Furthermore, it is well known in the machining process that, there is a positive correlation between the cutting edge recessing and the increasing of the cutting forces, highlighted for instance by [START_REF] Kivimaa | Cutting Force in Wood Working[END_REF] in the case of the woodworking, by [START_REF] Stewart | Tool forces and edge recession from cutting Medium Density Fiberboard[END_REF] in the case of MDF machining, and by [START_REF] Ravindra | Modelling of tool wear based on cutting forces in turning[END_REF] in the case of metal cutting. According to [START_REF] Sheikh-Ahmed | Measurement of tool wear and dulling in the machining of particleboard[END_REF], there is a close relationship between the value of the tangential forces divided by the normal forces and the tool lifespan. This gives the opportunity to evaluate the tool wear by measuring cutting forces as shown by [START_REF] Furet | Détermination de protocoles d'utilisation du capteur d'efforts en perc ¸age (Intellitool)[END_REF]. Furthermore, we can monitor the tool wear during the machining process and change the tool if necessary. To be aware of the cutting forces, according to [START_REF] Boucher | Influence of helix angle and density variation on the cutting force in wood-based products machining[END_REF], may also help in the field of tools design i.e. adapting the geometry of the tool, improve the tool holding during machining, and limit tool vibrations which could result in a better final wood surface quality.

The aim of this study is to development an on-line monitoring system dedicated to the routing of the Aleppo pine. By identifying one or more criteria to efficiently quantify and identify rapid tool wear in situ, it can be determined when to change the tool. Furthermore, the Aleppo pine wood is not enough valorised, it contains naturally many defects (knots, resins... etc.), more than other commercial pine wood (Douglas-fir, Norway spruce, Woodland pine...etc.). This makes it necessary to design new specific and strong tools. Our interest for Aleppo pine is that this wood species is widely introduced around Mediterranean Sea, and constitutes the first resource of wood in Algeria.

Materials and methods

Sampling

A CNC routing machine RECORD1 of SCM was utilised, fitted with a carbide tool from LEITZ with just a single cutting edge (z =1) (Fig. 1(a)). The cutting tool was designed to have two cutting edges, so to avoid out of balance effects caused by having only one insert, two inserts were fitted, one was deliberately blunted (so as not to cut) and the second was sharp. All cutting parameters were fixed (cutting geometry; feed rate per teeth...), they are listed in Table 1. The cutting direction related to the part anisotropy axis was 90 • -0 • (peripheral milling). Wood moisture content (MC) and temperature were controlled and fixed, respectively at 12% and 20 • C. The other parameters were selected in order to obtain prompt treatment and to identify the most appropriate method to analyze the forces. One single sample was prepared to dimensions of 24 mm × 24 mm × 3.5 mm and fixed on the 3 axes Dynamometer for the force acquisitions (Fig. 1(b)). In this way, the influence of wood variability on the cutting force acquisition was reduced. The tool wear was obtained by machining samples taken randomly from 19 Aleppo pine (Pinus halepensis) trees from the Aubagne region in France. These trees were 60 years old and representatives of the resource of this kind of wood around the Mediterranean Sea. Their average diameter was 36 cm (from 32 to 41 cm). The heart wood represent 90% of the wood used for this study, sap wood was removed by sawing in the first transformation. The mean density was: 0.55 g/cm 3 at 12% MC (from 0.53 to 0.58 g/cm 3 ). The other physical properties are summarized in

Cutting force measurement

A Kistler 9257A 3 axes Dynamometer was used to measure the cutting forces. The cutting forces were measured according to the two principal directions, tangential and normal, respectively FT and FN (Fig. 2). We define as the angle between cutting force F c and its tangential component FT. To reduce the effect of system vibrations, the data was filtered using a 4th order low pass filter set to 500 Hz. 500 Hz was selected because of the machining frequency being about 462 Hz. Data was computed using Dasylab software. All measured parameters are summarized in Table 3.

Wear measurement

Fig. 3 shows the observation direction and the wear criteria used in this study (wear mark Vb defined by Gottloeber et al. ( 2001)). The wear was observed using a LEICA MZ12 micro- scope, which had a suitable field depth (Fig. 4(a)). To minimize positional measurement errors, the observation of the wear was made directly on the cutting tool without removing the insert. A CCD camera, linked to the microscope, was used for the picture acquisition and LEICA software was used for the image computation and to enable on-line measurements. Wear measurements were made along the length of the active cutting edge and each value is an average of 10 measurements, an example of these measurements is shown in Fig. 4(b).

Results and discussion

The shape of the curve about the evolution of wear according to the cutting length is similar to the typical wear curve of the tools (Fig. 5). We can consider that the running period is between the beginning of the cutting process and 850 m of cutting length because of the relative fast wear. After 850 m, the tools wear is almost stable and the slope of the curve is less strong than in the first part. The observations are confirmed by the forces variations (Figs. 6 and7). The high forces variation observed for the 540 m cutting length point is due to the presence of "knots", the knots density is higher than the rest of the wood, which explains the increase in force values. This effect has been well reported by [START_REF] Eyma | Strains and cutting forces involved in the solid wood rotating cutting process[END_REF]. The knot is one of the natural wood defects; it is a branch that is included in the wood of a tree stem by growth around its base [START_REF] Kollmann | Principles of Wood Science and Technology[END_REF]. It is interesting to observe that the variation of the computed angle (defined earlier as the angle between cutting force F c and its tangential component FT) has also confirmed that up to 850 m cut length the cutting tool is still in the running period and after 850 m it can be assumed to pass to a stability period (Fig. 8). In fact, in the running zone is highly oscillating comparing to the stability zone where the oscillation tends towards zero.

Conclusions and prospects

This work provides some interesting results about the tool wear in a machining process of the North African Aleppo pine: the correlation between the tool wear and cutting forces shows that, the running period is about 850 m of cutting length. In this period, angle variation is unstable before it is around -1.1 • in the stability zone. This confirms the existence of the first two separate areas cited previously, that characterize the behaviour of the tool edge recession, which are running (abrupt wear) and the linear wear (known as stability period).

As a consequence we can say that angle is a good criterion for tool wear estimation in these practical machining conditions. We propose to use the same approach with other tools. Another important aspect, particularly in the case of wood, is to estimate the surface quality by measuring its roughness and wettability and to make the link with tool wear i.e. the transformation of wood by removal of matter generates new surfaces which have a precise functionality and will all receive protective films or glues. This highlights the need to characterise wood surface wettability. Depending on the various fields of application, a future aim of this study is to make practical recommendations to wood professionals, summarized in a good practice guide to enhance Aleppo pine solid woods value.

Fig. 1 -

 1 Fig. 1 -(a) Cutting Tool (80 mm diameter); (b) fixture of the sample on the Router table.

Fig. 2 -

 2 Fig. 2 -Distribution of the cutting forces.

Fig. 3 -

 3 Fig. 3 -Wear measurement principle.

Fig. 4 -

 4 Fig. 4 -(a) The system of observation of wear and (b) a cutting edge measurement example.

Fig. 5 -

 5 Fig. 5 -Wear evolution according to cutting length.

Fig. 8 -

 8 Fig. 8 -Evolution of the angle according to the cutting length.

Table 1 -

 1 Cutting parameters

	Designation	Values	Designation	Values
	Rake angle	10 •	Cutting depth	2 mm
	Sharpness angle 55 •	Rotational speed 6000 rpm
	Clearance angle 25 •	Cutting speed	25.12 m/s
	Tool material	98% WC + 2% Co Feed speed	5 m/min
	Tool diameter	80 mm	Moisture content 12%

Table 2 .

 2 T able2-Ph ysical properties of 19 Aleppo pine trees

	Designation	Values	Variations interval
	Diameter average	36 cm	From 32 to 41 cm
	Density average	0.55 g/cm 3 at 12%	From 0.53 to 0.58 g/cm 3
	Average Monnin	2.8	From 2.4 to 3.1
	hardness
	Young's modulus	11,200 MPa	8300 to 13,500 MPa
	average
	Bending fracture	82MP a	From63to98MP a
	stress

Table 3 -

 3 Cutting forces measurement and acquisition parameters

	Designation	Values	Designation	Values
	Dynamometer	3.5 kHz Filter frequency	500 Hz
	frequency
	Acquisition	10 kHz	Filter order	4th
	frequency
	Time of a single	0.01 s	Filter kind	Low pass
	rotation
	Tool frequency	100 Hz	Acquisition system 3 axes
				Dynamometer

Acknowledgements

This work could be performed because of the financial support of FCBA Institut Technologique, of FSE (Fond Social Européen), of CRB (Conseil Régional de Bourgogne) and of the University of Boumérdes.


B Darties

G Simonin

R Giroudeau

J.-C König

Coupled-tasks in presence of bipartite compatibilities graphs

We tackle the makespan minimization coupled-tasks problem in presence of incompatibility constraints. In particular, we focus on stretched coupled-tasks, i.e.coupled-tasks having the same sub-tasks execution time and idle time duration. We study several problems in the framework of classic complexity and approximation for which the compatibility graph is bipartite (star, chain, . . .). In such context, we design efficient polynomial-time approximation algorithms according to different parameters of the scheduling problem.

Introduction

We consider a non-preemptive coupled-tasks scheduling problem in presence of incompatibility constraint on a single processor. From the point of view of scheduling theory, the problem is also defined as a scheduling problem with exact delays on single machine. In this article, we will show the close relationship between coupled-task in presence of incompatibility constraint and the classic bin packing problem in the framework of complexity and approximation.

The coupled-tasks model, was first introduced by Shapiro [START_REF] Shapiro | Scheduling coupled tasks[END_REF] in order to model some data acquisition processes i.e. radar sensors: a sensor emits a radio pulse (first sub-task), and finally listen for an echo reply (second sub-task). Between these two instants (emission and reception), clearly there is an idle time due to the propagation, in both sides, of radio pulse. Therefore, a coupled-task is constituted by the triplet: the two sub-tasks and the idle between them. Thus, in order to minimise the makespan (schedule length), it is necessary to execute one or several different sub-tasks during the idle time of a coupled-task. Therefore, the aim is to find a best packing of coupled-tasks in which the sum of idle times is minimised. Notice that in the basic model, all coupled-tasks may be executed in each other according to processing time of sub-tasks and the duration of the idle time. Hereafter, we consider a relaxation of the previous model in which for a fixed coupled-task A there are only a subset of coupled-tasks compatible to be processed in the idle time of A. This model is motivated by the problem of data acquisition in presence of incompatibility constraint in a submarine torpedo. A collection of sensors acquires data for the torpedo. The incompatibility constraint is expressed to prevent interference issues caused by tasks using sensors working at the same frequency. So, the constraints are represented by a compatibility graph in which vertices are the coupled-tasks and edges represent compatibility between two tasks. In this article the variation of the complexity according to severals structural parameters are considered and some efficient polynomial-time approximation results on N P-hard instances are presented without omitting the relationship to bin packing problems.

Above all, we will show the close relationship between the studied problem and four packing-related problems, for which known approximation will be used as routine for scheduling coupled-tasks problem:

1. The subset sum (ss) problem: given a set S of n positive values and v ∈ IN, the aim is to find a subset S * ⊆ S such that i∈S * i = v. This problem is known to be N P-complete (see [START_REF] Garey | Computers and Intractability: A guide to the theory of NP-completeness[END_REF]). The optimization version is sometimes viewed as a knapsack problem, where each item profit and weight coincide to a value in S, the knapsack capacity is v, and the aim is to find the set of packable items with maximum profit. 2. The multiple subset sum (mss) problem: variant of bin packing in which a number of identical bins are given and one aims to maximize the overall weight of the items packed in the bins without violating the constraint on the capacity of each bin. The problem is a special case of the Multiple knapsack problem in which all knapsacks have the same capacity and the item profits and weights coincide. mss admits a PT AS [START_REF] Caprara | The Multiple Subset Sum Problem[END_REF] and a 3 4 -approximation algorithm [START_REF] Caprara | A 3/4-Approximation Algorithm for Multiple Subset Sum[END_REF], but does not admit a FPT AS even for only two knapsacks. 3. multiple subset sum with different knapsack capacities (mssdc) [START_REF] Caprara | A PTAS for the Multiple Subset Sum Problem with different knapsack capacities[END_REF] is an extension of mss considering different bin capacities. mssdc also admits a PT AS [START_REF] Caprara | A PTAS for the Multiple Subset Sum Problem with different knapsack capacities[END_REF]. 4. As a generalization of mssdc, multiple knapsack assignment restriction (mkar) problem consists in packing weighted items into non-identical capacity-constrained bins, with the additional constraint that each item can be packed into some bins only. Each item as a profit, the objective here is to maximize the sum of profits of packed items. Considering that the profit of each item is equal to its weight, [START_REF] Dawande | Approximation Algorithms for the Multiple Knapsack Problem with Assignment Restrictions[END_REF] proposed a 1 2 -approximation.

Presentation of coupled-tasks and related work

We model a task A i with a triplet (a i , L i , b i ), where a i (resp. b i ) is the duration of the first (resp. second) sub-task, and L i the idle time to respect between the execution of sub-tasks. We note A the set of tasks, and describe the incompatibility constraint between tasks with a graph G c = (A, E). There is an edge (A i , A j ) ∈ E iff a (or both) sub-task from A i may be scheduled during the idle time of A j or reciprocally. In a valid schedule, we said that A i is packed into A j if the entire task A i is scheduled during the idle time of A j . This is only possible when

a i + L i + b i ≤ L j . We call stretched coupled-task a task A i such that a i = L i = b i = α(A i ),
where α(A i ) is the stretch factor of task A i . And for any set W of tasks, we define seq(W ) = 3 x∈W α(x).

Due to the combinatorial nature of the parameters of the problem, we use the Graham's notation scheme α|β|γ [START_REF] Graham | Optimization and Approximation in Deterministic Sequencing and Scheduling: a Survey[END_REF] (respectively the machine environment, job characteristic and objective function) to characterize the problems related to coupled-tasks. The job characteristics summarizes the conditions made on the values of a i , L i , b i (independent between tasks, or equal to a constant), and the shape of the compatibility graph G. The coupled-tasks scheduling problems under incompatibility constraints has been studied in the framework of classic complexity and approximation in [START_REF] Simonin | Isomorphic coupled-task scheduling problem with compatibility constraints on a single processor[END_REF][START_REF] König | Complexity and approximation for scheduling problem for coupled-tasks in presence of compatibility tasks[END_REF].

3 Stretched coupled-task: model and contribution

Model

This paper focuses on stretched coupled-tasks. In the rest of the paper, all tasks are always stretched coupled-tasks, and, for two compatible tasks A j and A j to be scheduled in parallel, one of the following conditions must hold: 1. α(A i ) = α(A j ): the idle time of one task is fully exploited to schedule a sub-task from the other (i.e. b i is scheduled during L j , and a j is scheduled during L i ), and the completion of the two tasks is done without idle time. 2. 3α(A i ) ≤ α(A j ): task A i is fully executed during the idle time L j of A j .

From this observation, one can obtain from the compatibility graph G = (A, E) a directed compatibility graph G c = (A, E c ) by assigning a direction to each edge E from the task with the lowest stretch factor to the task with the highest one. If two compatible tasks x and y have the same stretch factor, then E c contains both the arc (x, y) and the arc inverted (y, x). Remark that if for any pair of compatible tasks x and y we have α(x) = α(y), then G c is a directed acyclic graph.

We note N G (v) the neighbourhood of v in G. We note d G (v) = |N (x)| the degree of v in G, and ∆ G the maximum degree of G. As we focus our work on bipartite graphs, we recall that a k-stage bipartite graph

is a digraph G = (V 0 ∪ • • • ∪ V k , E 1 ∪ • • • ∪ E k )
where V 0 . . . V k are disjoint vertex sets, and each arc in E i is from a vertex in V i to a vertex in V i+1 . The vertices of V i are said to be at rank i, and the subgraph

G i = (V i-1 ∪ V i , E i ) is called the i-th stage of G. For clarity, 1-stage bipartite graphs can be referred as triplet (X, Y, E) instead of (V 0 , V 1 , E).

Computational complexity

In this section,we present several N P-complete and polynomial results. We first show the problem is N P-hard even when the compatibility graph is a star (Theorem 1), but solvable with an O(n 3 ) time complexity algorithm when G is a chain (Theorem 2). Then we focus our analysis when G c is a 1-stage bipartite graph. We prove the problem is solvable with an O(n 3 ) polynomial algorithm if ∆ G = 2 (Theorem 3), but becomes N P-hard when ∆ G = 3 (Theorem 4).

Theorem 1. The problem 1|a

i = L i = b i = α(A i ), G = star|C max is
polynomial if the central node admits at least one outcoming arc.

-N P-hard if the central node admits only incoming arcs.

Proof

If there exists at least one outgoing arc (x, y) ∈ G c from the central node x, then the optimal solution consists in executing the x-task into the y-task, then in processing sequentially the remaining tasks after the completion of the y-task.

In the case where the central node admits only incoming arcs, first one can easily see that 1|α

(A i ) = a i = L i = b i , G = star|C max is N P. Second, we
propose the following polynomial construction from an instance of ss to an instance of our problem: ∀i ∈ S we add a coupled-task x with α(x) = i; let T be the set of these tasks; we add a task y with α y = a y = L y = b y = 3 × v; we define an incompatibility constraint between each task x ∈ T and y modelled by the compatibility graph G. In brief, G is a star with y as the central node. From this transformation, one can easily show the reduction between both problems.

Theorem 2. The problem 1|a

i = L i = b i = α(A i ), G = chain|C max admits a polynomial-time algorithm.
Sketch of proof Due to space limitation, we give only the main idea of the proof. Nevertheless, one can find the entire proof in the technical report [START_REF] Benoit Darties | Scheduling stretched coupled-tasks with compatibilities constraints : model, complexity and approximation results for some class of graphs[END_REF].

The proof consists first in simplifying the original instance by defining some elements of an optimal solution, in order to obtain a sub-instance where in any solution at most one task can be packed in another. This new instance of the problem can be solved in polynomial time by reducing it to the search of a minimum weighted perfect matching. Basically, this reduction consists in duplicating the compatibility graph G, then in linking each node to its clone with an edge. From this new graph, we add for each edge {x, y} a weight w({x, y}) corresponding to half the execution time of these two tasks, and to the processing time of x if y is the clone of x. Then we perform a minimum weighted perfect matching in O(n 2 m) by [START_REF] Edmonds | Maximum matching and a polyhedron with 0, 1 vertices[END_REF].

In following, we study the variation of the complexity in presence of a 1-stage bipartite graph according to the different values.

Theorem 3. The problem of deciding whether an instance of 1|a

i = L i = b i = α(A i ), G c =1-stage bipartite, ∆ Gc = 2|C max is polynomial.
Proof Let G c = (X, Y, E) be a 1-stage bipartite compatibility graph. Y -tasks will always be scheduled sequentially. The aim is to fill their idle time with a maximum of tasks of X, while the remained tasks will be executed after the Y -tasks. We just have to minimize the length of the remained tasks. Note that d Gc (y) ≤ 2. The algorithm use three steps :

1. for each task y ∈ Y such that 3 × α(x 1 ) + 3 × α(x 2 ) ≤ α(y) where x 1 and x 2 are the only two neighbors of Y , we add y to the schedule and execute x1 and x 2 sequentially during the idle time of y. Then we remove y, x 1 and x 2 from the instance. 2. Each remaining task y ∈ Y admits at most two incoming arcs (x 1 , y) and/or (x 2 , y). We add a weight α(x) to the arc (x, y) for each x ∈ N (y), then perform a maximum weight matching on G c in order to minimize the length of the remained tasks of X. Thus, the matched coupled-tasks are executed, and these tasks are removed from G c . 3. Then, remaining tasks from X are allotted sequentially after the other tasks.

The complexity of an algorithm is O(n 3 ).

Theorem 4. The problem of deciding whether an instance of 1|a

i = L i = b i = α(A i ), G c =1-stage bipartite, ∆ Gc = 3|C max
has a schedule of length at most 54n is N P-complete with n the number of tasks.

Proof It is easy to see that our problem is in N P. Our proof is based on a reduction from one-in-(2,3)sat(2, 1): does there exist an assignment of a set V of n boolean variables with n mod 3 ≡ 0, a set of n clauses of cardinality two and n/3 clauses of cardinality three such that:

-Each clause of cardinality 2 is equal to (x ∨ ȳ) for some x, y ∈ V with x = y.

-Each of the n literals x (resp. of the literals x) for x ∈ V belongs to one of the n clauses of cardinality 2, thus to only one of them. -Each of the n (positive) literals x belongs to one of the n/3 clauses of cardinality 3, thus to only one of them. -Whenever (x ∨ ȳ) is a clause of cardinality 2 for some x, y ∈ V, then x and y belong to different clauses of cardinality 3. We construct an instance π of our problem in following way (see Figure 1): 1. For all x ∈ V, we introduce four variable-tasks: x, x ′ , x and x′ with (a i , L i , b i ) = (1, 1, 1), ∀i ∈ {x, x ′ , x, x′ }. This variable-tasks set is noted VT . 2. For all x ∈ V, we introduce three literal-tasks L x , C x and Cx with L x = (2, 2, 2); C x = Cx = [START_REF] Edmonds | Maximum matching and a polyhedron with 0, 1 vertices[END_REF][START_REF] Edmonds | Maximum matching and a polyhedron with 0, 1 vertices[END_REF][START_REF] Edmonds | Maximum matching and a polyhedron with 0, 1 vertices[END_REF]. The set of literal-tasks is denoted LT .

3. For all clauses with a length of three, we introduce two clause-tasks C i and Ci with C i = (3, 3, 3) and Ci = (6, 6, 6). 4. For all clauses with a length of two, we introduce one clause-task C i with C i = (3, 3, 3). The set of clause-tasks is denoted CT . 5. The following arcs model the incompatibility constraints:

(a) For all boolean variables x ∈ V, we add the arcs (L x , C x ) and (L x , Cx ) (b) For all clauses with a length of three denoted C i = (y ∨ z ∨ t), we add the arcs (y, C i ), (z, C i ), (t, C i ) and (ȳ ′ , Ci ), (z ′ , Ci ), ( t′ , Ci ). (c) For all clauses with a length of two denoted C i = (x ∨ ȳ), we add the arcs (x ′ , C i ) and (ȳ, C i ). (d) Finally, we add the arcs (x, C x ), (x ′ , C x ) and (x, Cx ) and (x ′ , Cx ).

This transformation can be computed clearly in polynomial time. The proposed compatibility graph is 1-stage bipartite and d Gc (x) ≤ 3, ∀x ∈ VT ∪ LT ∪ CT .

In follows, we say that a task x is merged to a task y, if it exists a incompatibility constraint from x to y; i.e.the coupled-task x may be executed during the idle of coupled-task y.

⇒ Let us first assume that there is a schedule with length of 54n at most. We prove that there is a truth assignment I : V → {0, 1} such that each clause in C has exactly one true literal. We make some essentials remarks:

1. The length of the schedule is given by an execution time of the coupled-tasks admitting only incoming arcs, and the value is

54n = 3α CT |CT |+α LT (|LT |- |{L x , x ∈ V}|) = 9|{C i ∈ CT of length 2 and 3}| + 18|{ Ci ∈ CT }| + 18|{C x
and Cx ∈ LT }| = 9 × 4n 3 + 18 × n 3 + 18 × 2n. Thus, all tasks from VT ∪ {L x , x ∈ V} must be merged with tasks from CT ∪ (LT -{L x , x ∈ V}). 2. By the construction, at most three tasks can be merged together. 3. L x is merged with C x or Cx . 4. The allocation of coupled-tasks from CT ∪ (LT -{L x , x ∈ V}) leads to 18n idle time. The length of the variable-tasks VT and L x equals 18n (in these coupled-tasks there are 6n idle times). 5. If the variable-tasks x and x ′ are not merged simultaneously with C x , i.e.only one of these tasks is merged with C x , so, by with the previous discussion, it is necessary to merge a literal-task L y , with x = y one variable-task (ȳ or ȳ′ ) with C y or Cy . It is impossible by size of coupled-tasks. In the same ways, the variable-tasks x et x′ are merged simultaneously with Cx . 6. Hence, first x and x ′ are merged with C x or with clause-task where the variable x occurs. Second, x and x′ are merged with Cx or a clause-task.

So, we affect the value "true" to the variable l iff the variable-task l is merged with clause-task(s) corresponding to the clause where the variable l occurs. It is obvious to see that in the clause of length three and two we have one and only one literal equal to "true".

⇐ Conversely, we suppose that there is a truth assignment I : V → {0, 1}, such that each clause in C has exactly one true literal.

1. If the variable x = true then we merged the vertices L x with C x ; x with the clause-task C i corresponding to the clause of length three which x occurs;

x ′ with the clause-task C i corresponding to the clause of length two which x occurs; and x, x′ with Cx . 2. If the variable x = f alse then we merged the vertices L x with Cx ; x with the clause-task corresponding to the clause of length two which x occurs; x′ with the clause-task Ci corresponding to the clause (C) of length three which x occurs; and x, x ′ with C x .

For a feasible schedule, it is sufficient to merge vertices which are in the same partition. Thus, the length of the schedule is at most 54n.

5 Polynomial-time approximation algorithms

Star graph

Theorem 5.

1|a i = L i = b i = α(A i ), G = star|C max admits a FPT AS.
Proof We may use the solution given by the subset sum (ss) (see [START_REF] Ibarra | Fast approximation algorithms for the Knapsack and Sum of Subset problems[END_REF] and [START_REF] Kellerer | An efficient fully polynomial approximation scheme for the subset-sum problem[END_REF])). Indeed, the schedule is follows: first the central node is executed, second during its idle time we process the coupled-tasks chosen by an FPT AS algorithm from ss, and finally the remaining tasks are processed after the completion of the central node.

1-stage bipartite graph

Scheduling coupled-tasks during the idle time of others tasks can be related to packing problems, especially when the compatibility graph G c is a bipartite graph. In the following, we propose several approximation results when G c is a 1-stage bipartite graph.

Lemma 1. Let P be a problem with P ∈ {mkar mssdc, mss} such that P admits a ρ-approximation, then the following problems

1. 1|a i = L i = b i = α(A i ), G c = 1-stage bipartite|C max , 2. 1|α i = a i = L i = b i , complete 1-stage bipartite|C max 3. 1|α i = a i = L i = b i , complete 1-stage bipartite|C max
where the compatibility graph is a complete bipartite G=(X,Y), and all the tasks from Y have the same α(y). are approximable to a factor 1 + (1-ρ) 3 . E)and a stretch factor function α : X ∪ Y → IN. In such instance, any valid schedule consists in finding for each task y ∈ Y a subset of compatible tasks X y ⊆ X to pack into y ∈ Y , each task of x being packed at most once. Let X p = ∪ y∈Y X y be the union of tasks of X packed into a task from Y , and let X p be the set of remaining tasks, with X p = X/X p . Obviously, we have:

Proof 1. Let consider an instance of 1|α

i = a i = L i = b i , G c =1-stage bipartite|C max with G c = (X, Y,
seq(X p ) + seq(X p) = seq(X) (1)
As Y is an independent set in G c , tasks from Y have to be scheduled sequentially in any (optimal) solution. The length of any schedule S is then the processing time of Y -tasks plus the execution time of the X p-tasks. Formally:

C max (S) = seq(Y ) + seq(X p) = seq(Y ) + seq(X) -seq(X p ). (2)
We use here a reduction to mkar: each task x from X is an item having a weight 3.α(x), each task from Y is a bin with capacity α(y), and each item x can be packed on y if and only if the edge {x, y} belongs to G c . Using algorithms and results from the literature, one can compute the set X p of packed items. The cost of the solution for the mkar problem is seq(X p ).

If mkar is approximable to a factor ρ, then we have:

seq(X p ) ≥ ρ × seq(X * p ), (3)
where X * p is the set of packable items with the maximum profit. Combining Eq. ( 2) and ( 3), we obtain a schedule S with a length equal to:

C max (S) ≤ seq(Y ) + seq(X) -ρ × seq(X * p ) (4)
As X and Y are two fixed sets, an optimal solution S * with minimal length C max (S * ) is obtained when seq(X p ) is maximum, i.e. when X p = X * p . Therefore, the ratio obtained between our solution S and the optimal one S * is:

C max (S) C max (S * ) ≤ seq(Y ) + seq(X) -ρ × seq(X * p ) seq(Y ) + seq(X) -seq(X * p ) ≤ 1+ (1 -ρ) × seq(X * p ) seq(Y ) + seq(X) -seq(X * p ) (5)
By definition, X * p ⊆ X. Moreover, as the processing time of X * p cannot excess the idle time of tasks from Y , we obtain: seq(X * p ) ≤ 1 3 seq(Y ). And thus combined to Eq. ( 5), we obtain the desired upper bound:

C max (S) C max (S * ) ≤ 1 + (1 -ρ) 3 (6)
2. For the problem 1|α i = a i = L i = b i , complete 1-stage bipartite|C max , the proof is identical using mssdc as a special case of mkar where each item can be packed in any bin.

For the problem 1|α

i = a i = L i = b i , complete 1-stage bipartite|C max
where all the tasks from Y have the same stretch factor α(y), the proof is identical as previously since mssdc is a generalisation of mss.

Theorem 6. These problems admit a polynomial-time approximation algorithm:

1. The problem 1|a i = L i = b i = α(A i ), G c = 1-stage bipartite|C max is ap- proximable to a factor 7 6 . 2. The problem 1|a i = L i = b i = α(A i ), G c = complete 1-stage bipartite|C max admits a PT AS. 3. The problem 1|a i = L i = b i = α(A i ), G c = complete 1-stage bipartite|C max ,
where all the tasks from Y have the same stretch factor α(y):

(a) is approximable to a factor 13 12 . (b) admits a PT AS. Proof 1. Authors from [START_REF] Dawande | Approximation Algorithms for the Multiple Knapsack Problem with Assignment Restrictions[END_REF] proposed a ρ = 1 2 -approximation algorithm for mkar. Reusing this result with Lemma 1, we obtain a 7 6 -approximation. 2. We know that mssdc admits a PT AS [START_REF] Caprara | A PTAS for the Multiple Subset Sum Problem with different knapsack capacities[END_REF], i.e. ρ = 1-ǫ. Using this algorithm to compute such a PT AS and the Lemma 1, we obtain an approximation ratio of 1 + ǫ 3 for this problem. 3. (a) Authors from [START_REF] Caprara | A 3/4-Approximation Algorithm for Multiple Subset Sum[END_REF] proposed a ρ = 3 4 -approximation algorithm for mss. Reusing this result and the Lemma 1, we obtain a 13 12 -approximation. (b) They also proved that mss admits a PT AS [START_REF] Caprara | The Multiple Subset Sum Problem[END_REF] , i.e. ρ = 1ǫ. Using the algorithm to compute such a PT AS and the Lemma 1, we obtain an approximation ratio of 1 + ǫ 3 .

2-stage bipartite graph

Theorem 7. The problem 1|a i = L i = b i = α(A i ), G c = 2-stage bipartite|C max is approximable to a factor 13 9 .

Proof We consider an instance of the problem whith

G c = (V 0 ∪V 1 ∪V 2 , E 1 ∪E 2 )
, where each arc in E i is oriented from a vertex in V i to another one in V i+1 , for i ∈ 1, 2. Before presenting our heuristic and the analyse of its approximation factor, we will give several notations, properties and equations in relation with the specificities of this instance, in any (optimal) solution:

-∀i = 0, 1, let V ip (p=packed), (resp. V ia (a=alone) ) be the set of tasks merged (resp. remaining) into any task from V i+1 in a solution S, and V ib (b=box) the set of tasks scheduled with some tasks from V i-1 merged into it. This notation is extended to an optimal solution S * by adding a star in the involved variables. -Given any solution S to the problem and considering the specificities of the instance, note that

{V 0p , V 0a } is a partition of V 0 , G c , {V 1p , V 1a , V 1b } is a partition of V 1 , and G c , {V 2a , V 2b } is a partition of V 2 .
-Any solution would consists in scheduling first each task with at least one task merged into it, then to schedule the remaining tasks (alone). Given an optimal solution S * , the length of S * is given by the following equation:

S * = seq(V 1 * b ) + seq(V 2 b) + seq(V 0 * a ) + seq(V 1 * a ) + seq(V 2 * a ) S * = seq(V 2 ) + seq(V 1 * b ) + seq(V 0 * a ) + seq(V 1 * a ) (7)
One can remark that V 0 * p and V 1 * p are not part of the equation, as they are scheduled during the idle time of

V 1 * b and V 2 * b . -Let consider an restricted instance of G c to a sub-graph G 0 = G c [V 0 ∪ V 1 ] (resp. G 1 = G c [V 1 ∪ V 2 ]) which is the 1-th (resp. 2-th) stage of G c . Let S[G 0 ] (resp. S * [G 0 ]) be any (an optimal) solution on G 0 , V 0p [G 0 ] (resp. V 0 * p [G 0 ]) is the set of tasks from V 0 packed into tasks from V 1 in S[G 0 ] (resp. S * [G 0 ]), and V 0a [G 0 ] (resp. V 0 * a [G 0 ]
) the set of remaining tasks. In addition tho these notation, let V 1b [G 0 ] be the set of tasks from V 1 with at least one task from V 0 merged into them, and V 1a [G 0 ] the remaining tasks. A first observation gives for G 0 :

S

* [G 0 ] = seq(V 1 ) + V 0 * a [G 0 ] (8)
-From Theorem 6, Lemma 1, and the demonstration presented in their proof from [START_REF] Dawande | Approximation Algorithms for the Multiple Knapsack Problem with Assignment Restrictions[END_REF], several equations can be computed for a solution S[G 0 ]:

seq(V 0p [G 0 ]) ≥ 1 2 seq(V 0 * p [G 0 ]) (9)
seq(V 0 * a [G 1 ]) ≤ seq(V 0 * a ) (10)
seq(V 0p [G 0 ]) + seq(V 0a [G 0 ]) = seq(V 0 * p [G 0 ]) + seq(V 0 * a [G 0 ]) = seq(V 0 ) (11) seq(V 0a [G 0 ]) ≤ seq(V 0 * a [G 0 ]) + 1 2 seq(V 0 * p [G 0 ]) ≤ seq(V 0 * a ) + 1 2 seq(V 0 * p [G 0 ]) (12)
-We use an analog reasoning on the sub-graph G 1 with equivalent notations for V 1 and V 2 , and we obtain:

seq(V 1p [G 1 ]) ≥ 1 2 seq(V 1 * p [G 1 ]) (13)
seq(V 1a [G 1 ]) ≤ seq(V 1 * a [G 1 ])+1/2seq(V 1 * p [G 1 ]) ≤ seq(V 1 * a )+1/2seq(V 1 * p [G 1 ]) (14)
From this notations and observation, we can propose a good heuristic. We design the feasible solution S for G c as follows:

-We compute a solution S[G 1 ] on G 1 , then we add to S each task from V 2 and the tasks from

V 1 merged into them (i.e. V 1p [G 1 ]) in S[G 1 ]. -Then we compute a solution S[G 0 ] on G 0 , then we add to S each task v from V 1b [G 0 ]/V 1p [G 1 ] and the tasks from V 0 merged into them. -Tasks V 1a [G 1 ]/V 1b [G 0 ] and V 0a [G 0 ] are added to S sequentially.
-We note V conf lict the set of remaining tasks, i.e. the set of tasks from V 0 which are merged into a task

v ∈ V 1 in S[G 0 ], thus that v is merged into a task from V 2 in S[G 1 ].
Remark that:

seq(V 1b [G 0 ]/V 1p [G 1 ]) + seq(V 1a [G 1 ]/V 1b [G 0 ]) = V 1a [G 1 ] (15)
Thus the cost of our solution S is

S = seq(V 2 ) + seq(V 1a [G 1 ]) + seq(V 0a [G 0 ]) + seq(V conf lict ) (16)
It is also clear that:

seq(V conf lict ) ≤ 1 3 seq(V 1p [G 1 ]) ≤ 1 3 seq(V 1 * p [G 1 ]) (17)
Using Equations ( 12), ( 14) and (17) in Equation ( 16), we obtain

S ≤ seq(V 2 ) + seq(V 1 * a ) + 5 6 seq(V 1 * p [G 1 ]) + seq(V 0 * a ) + 1 2 seq(V 0 * p [G 0 ]) (18) ≤ S * + 5 6 seq(V 1 * p [G 1 ]) + 1 2 seq(V 0 * p [G 0 ]) , using Equation (7) (19)
We know that S * ≥ seq(V 2 ) and S * ≥ seq(V 1 ), as V 1 is an independent set of G c . We also know that tasks from

(V 1 * p [G 1 ]) (resp. (V 0 * p [G 0 ]
)) must be merged into tasks from V 2 (resp. V 1 ) and cannot exceed the idle time of V 2 (resp. V 1 ), implying that seq(

V 1 * p [G 1 ])) ≤ 1 3 seq(V 2 ) (resp. seq(V 0 * p [G 0 ])) ≤ 1 3 seq(V 1 )
). One can write the following :

5 6 seq(V 1 * p [G 1 ]) S * ≤ 5 6 × 1 3 seq(V 2 ) seq(V 2 ) ≤ 5 18 (20)

Conclusion

The results proposed in this paper are summarised in Table 1. New presented results suggest the main problem of coupled tasks scheduling remains difficult even for restrictive instances, here stretched coupled-tasks when the constraint graph is a bipartite graph. When we consider stretched coupled-tasks, the maximum degree ∆ G seems to play an important role on the problem complexity, as the problem is already N P-Hard to solve when the constraint graph is a star. Approximation results presented in this paper show the problem can be approximated with interesting constant ratio on k-stage bipartite graphs for k = 1 or 2. The presented approach suggests a generalisation is possible for k ≥ 3. This part constitutes one perspective of this work. Other perspective would consists to study coupled-tasks on other significant topologies, including degree-bounded trees, or regular topologies like the grid.

Fig. 1 .

 1 Fig. 1. A partial compatibility graph for the N P-completeness of the scheduling problem 1|bipartite of depth one, d(Gc) ≤ 3, αi = ai = Li = bi|Cmax

  Equations (19), (20) and (21) the proof is concluded:

Acknowledgment

This work has been funded by the regional council of Burgundy.

Topology

Complexity

Approximation uug(Gc)=Star graph N P -C (Theorem 1) F PT AS (Theorem 5) uug(Gc)=Chain graph

N P -C (Theorem 4) 7 6 -APX (Theorem 6) Gc= complete 1-stage bipartite N P -C (see [START_REF] König | Complexity and approximation for scheduling problem for coupled-tasks in presence of compatibility tasks[END_REF]) PT AS (Theorem 6) Gc= complete 1-stage bipartite N P -C (see [START_REF] König | Complexity and approximation for scheduling problem for coupled-tasks in presence of compatibility tasks[END_REF]) PT AS (Theorem 6) with constraint α(x) = α(y), ∀x, y ∈ X1


Chen-Iafrate-Savin Model

Julien Toulouse

Andreas Savin

Carlo Adamo
email: adamo@ext.jussieu.fr

Validation and assessment of an accurate approach to the correlation problem in density functional theory: The Krieger

In the present paper, we validate and assess a meta-GGA correlation functional (KCIS), whose form and parameters are entirely derived only from first-principles criteria. In particular, we have carried out a detailed comparison with the most common, parameterized correlation functionals. Next, we propose a new model in which the KCIS correlation is integrated in a hybrid Hartree-Fock/Density Functional Theory scheme. In such approach only one, or two in the G2-optimized version, parameters are adjusted on experimental data, all the others being derived from purely theoretical considerations. The results obtained for a set of molecular properties, including H-bonded complexes, proton transfer model, S N 2 reaction and magnetic (EPR) properties, are satisfactory and comparable, if not better, with those delivered by the most common functionals including heavy parameterization. The way in which the whole functional is derived and the few empirical parameters used make the new exchange-correlation functional widely applicable.

Introduction

Density Functional Theory (DFT), combining good performances and low computational costs, is nowadays one of the most powerful and reliable tools of quantum chemistry for the computation of the electronic structure of atoms, molecules or solids [START_REF] Kock | A Chemist's Guide to Density Functional Theory[END_REF]. In the application of this theory within the Kohn-Sham formalism, the total electronic energy can be exactly expressed as a functional of the electron density, the only contribution that needs to be approximated being the exchange-correlation term [START_REF] Parr | Density Functional Theory of Atoms and Molecules[END_REF]. Finding improved approximations to this contribution is therefore of vital importance in chemical applications to obtain accurate numerical molecular properties [START_REF] Burke | Electronic density functional theory: recent progress and new derivations[END_REF][START_REF] Adamo | [END_REF].

A large number of approximations have been developed, ranging from the simple local density approximation (LDA) involving only the electron density  to the generalized gradient approximations (GGA's) where the gradient of the density   has been introduced

to model the inhomogeneous character of molecular and real chemical systems. More recently, a new generation of functionals, called meta-GGA's and including the laplacian of the density (  2  ) and/or the kinetic energy density (  ), have attracted increasing attention.

Introducing more semi-local or non-local information compared to previous approximations, these meta-GGA functionals are very promising but have not been fully explored yet [5].

Two extreme strategies can be defined in designing new approximations. The first approach consists in proposing an empirical formula containing some parameters fitted on experimental data, with the related danger of obtaining functionals adapted only to the applications for which they have been developed. The second way is represented by designing a mathematical form which satisfies some (eventually the largest possible number of) rigorous theoretical conditions of the exact exchange-correlation functional. This approach does not assure the "chemical" applicability of the developed functional. Of course the two philosophies are compatible and can be mixed, as recently showed [6,7].

In practice, any functional is arbitrarily separated into an exchange and a correlation contributions, and approximated forms are often developed in an independent way. For instance, a large number of approximations have been proposed for the exchange term representing the major part of the combined contribution to the total electronic energy (see for instance reference [8][9][10][11][12]. In contrast, relatively few forms have been proposed for correlation functionals, perhaps due to the greater difficulty to represent this contribution. Here we recall some of the most common correlation forms, like those developed by Perdew and co-workers (like the P86, PW91 or PBE) [13][START_REF]PW91 correlation[END_REF][START_REF] Perdew | [END_REF], which are based on a GGA philosophy, while, only more recently, meta-GGA correlations appeared [16][17][18][19] Aside from these expressions, all rooted on the free-electron gas, several functionals based on the Colle-Salvetti approach [20] have been proposed, where the inhomogeneity of the systems is introduced via the second-order reduced density matrix [21][22][23]. The most known derivation has been developed by Lee, Yang and Parr and it is one of the most common and successful correlation functionals (usually referred to as LYP) [21].

In this context, Krieger, Chen, Iafrate and Savin have recently proposed a new meta-GGA approach (hereafter referred to as KCIS), based on the idea of a uniform electron gas with a gap in the excitation spectrum [START_REF] Krieger | Electron Correlations and Materials Properties[END_REF][START_REF] Rey | [END_REF]. The KCIS functional is, in our opinion, particularly appealing since it is built only from physical considerations, it preserves many of the know properties of the exact correlation energy (such as the so-called PBE conditions [START_REF] Perdew | [END_REF]) and it has no empirical parameters. Despite its promising features, the KCIS functional has never been tested beyond atomic or small molecular systems [START_REF] Krieger | Density functional theory and its applications to materials[END_REF].

The aim of this paper is, therefore, to explore the limits of the applicability of such a functional to chemical systems and to compare the KCIS performances with those provided by the most common correlation functionals. To this end, several exchange functional have been chosen as counter part either in pure DFT approaches or in a hybrid HF/DFT schemes [START_REF] Becke | [END_REF]. As in our previous works [7,10], we have tested the exchange-correlation functionals obtained in these ways by computing properties concerning a number of chemical systems or processes involving both covalent and weak interactions.

Exchange and Correlation Functionals used.

The KCIS correlation functional has the following expression [5]:

    r d F E c                     , , , , , , (1) where
) , 0 , ( ) , , ( ,                           GGAGAP c W GGAGAP c F (2)
In equation (2), GGAGAP c  is the correlation energy per particle of a uniform electron gas with a gap in the excitation spectrum and including gradient corrections, and W   is the Weizsäcker kinetic energy density:

) ( ) ( 8 1 ) ( 2 r r r         W (3)
Correlation functionals usually transfer the correlation energy from the uniform electron gas to the real system, eventually adding some gradient corrections. However, these two systems are very different: the spectrum of the uniform electron gas is a continuum of states which are occupied up to the Fermi level and unoccupied above. On the contrary, the molecular or atomic spectrum is made of discrete states which are occupied up to the HOMO orbital and unoccupied up to the ionization threshold, followed by a continuum of unoccupied states. It is therefore more judicious to model the correlation energy by starting from a uniform electron gas with a gap in the excitation spectrum roughly corresponding to the difference between the HOMO energy and the ionization threshold of the real system. It has been chosen to express this gap as the simplest semi-local functional positive everywhere and reducing to the ionization potential in the large r limit:

  2 2 8 1      G (4)
Furthermore, GGAGAP c  has been built so as to the three PBE theoretical conditions are fulfilled by the KCIS functional. Besides, the second term in expression (2) ensures the selfinteraction correction (SIC) of the functional which enables the correlation energy to vanish in the case of one-electron systems [28].

Since, any form for exchange functional have been proposed for the KCIS correlation, we have tested in this work several combinations using already existing exchange functionals. It is worth to notice that these tests are essential because the arbitrary association between an exchange and a correlation functional, developed separately, can in principle lead to any kind of numerical results. In fact, the quality of exchange-correlation functionals is almost always due to quite unpredictable compensations of errors of each contribution. Among the most common GGA approaches, we have chosen the Perdew-Wang 91 (PW91) functional [29], the Perdew-Burke-Ernzerhof (PBE) form [START_REF] Perdew | [END_REF] and the widely-used Becke 88 (B) exchange [9]. Concerning the meta-GGA approach, only the exchange part of the Perdew-Kurth-Zupan-Blaha (PKZB) model was tested [18]. Furthermore, the KCIS correlation functional has been also cast in a hybrid model of the general formula [START_REF]Becke one parameter hybrid[END_REF]:

  BKCIS xc B x exact x hyb xc E E E a E    (5)
The constant a can be fixed a priori to 1/4 by using theoretical considerations, or fitted to experimental data [START_REF]Becke one parameter hybrid[END_REF][START_REF] Adamo | [END_REF][32]. When using the Becke 88 exchange, the former approach using the theoretically determined value will be named B0KCIS, whereas, when using a fitted value for a, will be referred to as B1KCIS.

Computational details

All the DFT computations have been carried out within the unrestricted Kohn-Sham formalism as implemented in our development version of the Gaussian 99 package [START_REF] Frisch | Gaussian 99 Development version[END_REF]. We have added the KCIS correlation functional, as well as its first and second derivatives with respect to all the required variables so that all the available exchange functionals as well as all the standard features of the package (as geometry optimization, computation of excitation energies, magnetic properties) can be used in conjunction with KCIS.

A number of different basis sets have been used [START_REF] Frisch | Gaussian 98 User's Reference[END_REF]. In particular, the 6-311G(d,p) basis set of Pople and co-workers has been used to optimize all the molecular structures, since previous experience showed that a polarized valence triple- basis set generally provides nearly converged structural parameters by DFT methods. An extended basis set, namely the 6-311++G(3df,3pd) one, has next been used to evaluate all the energetic parameters (atomization and dissociation energies). Finally, magnetic properties have been evaluated

with the EPR-III basis set which has been purposely tailored for this kind of DFT applications.

When necessary, some computations have been also carried out with standard DFT model, such B3LYP, PBE or PBE0 [START_REF] Frisch | Gaussian 98 User's Reference[END_REF].

Results and discussion

The validation of a functional is still a problematic ground, where some arbitrary (and limited) choices have to be done. As matter of fact, that it is not possible a priori to predict the reliability of a chosen functional form in chemical applications and an investigation is then compulsory before routine applications. For instance, the B3LYP functional can be considered as an enlightening example for the large number of papers assessing its validity (see for instance the test cases reported in reference 1). Even if a large part of functionals have been tested only on thermochemistry (using the standard G2 molecular data set), some of them have been the subjects of more detailed investigations on large number of already wellknown chemical systems or processes [7,10,11]. Following this latter viewpoint, we have tested the quality of results given by the new functionals including KCIS, concerning thermodynamics, geometrical, spectroscopic and magnetic properties on chemical systems including some traditional pathological cases for DFT methods.

Atomic correlation energies

The first step for assessing the performances of the KCIS functional is to compared the correlation energies of the atoms belonging to the first two rows of the periodic table [START_REF] Chakravorty | [END_REF]. As can been seen by the results, obtained using HF densities, reported in table I, the accuracy of KCIS is at the same level than those of the most accurate correlation functionals, such as LYP, PKZB and PW91. This is a remarkable result, since these last functionals all contain parameters adjusted to experimental data, while this is not the case for the KCIS functional.

A part for this overall behaviour, two points must be stressed. First of all, it is noteworthy that the larger errors are obtained with atoms having an occupation of p orbitals. Indeed, the correlation energy of this kind of electronic configurations often involved contributions from virtual states which should be localized in the gap of the excitation spectrum but which are totally neglected by KCIS [XX]. As second point, we notice that although all the functional give similar errors, the mean absolute errors (mae's) are different for the first and second rows, the KCIS and the LYP functionals providing a more balance description between first and second rows atoms.

Atomization energies

In order to carry out a detailed selection of the numerous models that can be build with KCIS, we have computed for each functional the mae's with respect to experimental values on atomization energies of the molecules belonging to the standard G2-1 (55 molecules) and the extended G2 (148 molecules) sets [36][37][38]. These sets can be considered as a DFT "classical" test, representatives of covalently bonded systems [39]. Compared to the G2-1 set which contains small saturated molecules, the G2 set includes bigger molecules, eventually nonsaturated.

The results collected in table II, show that, among the most common functionals, the reference is the hybrid B3LYP model, with a mae of 2.4 kcal/mol over the G2-1 set and 3.1 kcal/mol for G2. In contrast, the PBE0 functional, which has no empirical parameters, leads to a relatively small mae for this kind of systems with 3.1 kcal/mol for the G2-1 set. The two most accurate model including the KCIS correlation are the PKZBKCIS and BKCIS with mae's of 4.1 and 8.2 kcal/mol over the G2-1 and G2 sets, respectively. The optimization of the parameter present in the exchange part of PKZB (D = 0.113) [18] leads to the PKZB opt KCIS functional (D opt = 0.104) which is not significantly better than the original PKZBKCIS (mae of 3.9 kcal/mol over the G2-1 set). Furthermore, the introduction of 25 % of exact exchange (PKZB0KCIS) only deteriorates the quality of the results (mae of 6.9

kcal/mol for G2-1). The construction of a more accurate model rather requires hybrid models with the BKCIS functional. In fact, upon introduction of 25 % of HF exchange into this functional (leading to the B0KCIS model), the mae falls to 3.0 kcal/mol for G2-1. The optimization of the percentage of exact exchange over the G2-1 set (B1KCIS with 23.9 % of exact exchange) does not induce significant improvements (see last line in table II). These good performances of B0KCIS, whereas the B0LYP functional is significantly less accurate, are, of course, to be credited to the different correlation.

As last test, we have carried out some computations on the G2 set, without the SCF procedure, using HF or PBE densities and only B0KCIS energy evaluation. This test is particularly interesting, since some recent implementations of exchange-correlation functionals do not include the derivatives need to the SCF procedure [6,18]. Our results (not reported in table II) show that with a PBE density the mae is 9.1 kcal/mol and it rises up to 11.7 kcal/mol with the HF densities. Both results are really far from the 3.0 kcal/mol obtained with the complete SCF procedure, thus warning about the use of auto-coherent procedure.

In the light of these first results, we have decided to study in details the performances of the KCIS functional, and especially those provided by the B0KICS model, which seems to be very accurate for covalent systems, even if it includes only one empirical parameter (in the exchange contribution). The natural reference to compare this model is, in our opinion, the similar B0LYP functional, but, in order to have a more complete picture, we consider also the popular B3LYP and PKZB functionals.

Atomic total energies and geometry optimization

As a first step, we have evaluated the total energies of some atoms, H through Ar. The results are collected in table III and are compared with accurate post-HF values (labeled as "exact" in the table) [START_REF] Frisch | Gaussian 99 Development version[END_REF]. It is interesting to evaluate the accuracy of functionals concerning atomic total energies since this property has been recently added to a wide training set used to test and optimize DFT models [40]. Over these 18 atoms, B0KCIS provides a mae of 0.008 Hartree.

This very good result are equaled only by B0LYP containing though four extra empirical parameters in the correlation contribution. We can also notice the astonishing accuracy of B0KCIS over the second row of the periodic table (mae of 0.004 Hartree).

After atomization and atomic energies, we have tested the ability of our new model to reproduce correct molecular geometrical parameters. Figure 1 reports the mae's over bond lengths of 32 molecules belonging to the G2 set computed with several functionals. With a mae of 0.008 Å, B0KCIS leads to about the same accuracy of other hybrid models PBE0, B3LYP and B0LYP. In fact, it is a well-known behavior: the introduction of a fraction of exact exchange in a functional generally leads to very accurate models concerning bond lengths. Indeed, since bond lengths are usually too short at the Hartree-Fock level and to long for GGA functionals, hybrid models benefit from a compensation of errors [START_REF] Adamo | [END_REF].

Other tests

Even if the exchange-correlation term to be approximated in the KS approach is usually decomposed in well-separated exchange and correlation parts, it is difficult to determine the dependence of one molecular properties on just one of these contributions. Few exceptions are represented by the van der Waals interactions and excitations dimers, both essentially ruled by the asymptotic behavior of the exchange functional [10,41]. Here, instead, we focused on more correlation-dependent properties/systems. In this context, small H-bonded dimers are an invaluable test, both for their chemical role and for the difficult to describe the subtle physical effects (electrostatic, charge transfer, polarization) present. It is therefore interesting to estimate the quality of DFT functionals on the kind of systems. Table VI collects geometrical parameters and interaction energies of three small H-bond dimers, namely (H 2 O) 2 , (HCl) 2 and (HF) 2 . The computed interaction energies have been corrected for BSSE error which has an order of magnitude of about 0.2 kcal/mol for the first two dimers with the large 6-311++G(3df,3pd) basis set, while this error is even lower for (HF) 2 . For the three systems, only the global energy minimum, of C s symmetry, has been investigated, even if other stable rearrangements could actually exist [42].

The water dimer has been studied in considerable detail and the currently accepted values for the binding energy and the OO distance are 5.4 ± 0.7 kcal/mol (not including ZPE corrections) and 2.952 Å, respectively [43,44]. Refined post-HF computations suggest that the interaction energy is close to 4.7 kcal/mol and give a shorter intermolecular distance (2.925 Å) [42,45]. As already well documented (see reference 1 for an almost complete bibliography), the equilibrium geometry of the water dimer is quite accurately reproduced by standard DFT methods, all the predicted OO distances falling between 2.87 and 2.91 Å, i.e.

close to the MP2 estimate [START_REF] Kock | A Chemist's Guide to Density Functional Theory[END_REF]. Our results are reported in table IV, and compared with experimental and post-HF data [42][43][44]. The most striking feature of these results is the B0KCIS distance, which is slightly longer (2.91 Å) than other DFT values, but closer to the experimental value (2.95 Å) [42]. At the same time the binding energy is slightly underestimated (3.8 vs. 4.7 kcal/mol for the best ab initio value). Interesting both B3LYP and B0LYP provide higher interactions energies and shorter intermolecular distances, thus

showing the correlation between these two quantities.

Similar trends are present for all the others two dimers, the B0KCIS providing always lower interactions energies and higher H-bridge distance than the B3LYP and B0LYP functionals.

In the case of the hydrogen fluoride dimer, the combined effects of shorter distance/lower energy results in a better agreement with the experimental data (see table IV) [47]. This is not the case for (HCl) 2 , where all the methods provide an poor accord with experiments [48]. A detailed comparison with literature data, suggests that more subtle effects, involving either exchange and correlation contributions, are at work in such H-bonded dimers [7].

The next test is represented by bimolecular nucleophilic substitution reactions (S N 2), which have an important place in organic chemistry. The symmetric Walden inversion, in particular, a excellent well-characterized example [START_REF] Shaik | Theoretical aspects of physical organic chemistry. The SN 2 mechanism[END_REF][START_REF] Truong | [END_REF][51][52][53][54]:

Cl -+ CH 3 Cl ClCH 3 +Cl -
This reaction has a double-well energetic profile, with two minima corresponding to the formation of a pre-and a post-reaction ion-molecule complex and a saddle point (SP) with

D 3h symmetry ([Cl … CH 3 … Cl] -) .
The most significant thermodynamic quantities are the complexation energy of the ion-molecule complex (E comp ), the activation energy, i.e. the relative energy of the D 3h saddle point with respect to the ion-molecule complex (E # ), and the overall barrier (E ovr ), defined as the difference between these two energies. The computed energies are collected in Table V. Apart from PBE0 which leads to a realistic reaction barrier, it is striking to notice that the more popular semi-empirical approaches, B3LYP or B0LYP, as well as more theoretically developed functionals, like PBE or PKZB, lead to zero or even negative values for E ovr . Whereas the computation of initial closed-shell reagents does not generate particular difficulties, the determination of the energy of the charged transition state energy [Cl … CH 3

…

Cl] -by DFT approaches is more involved [52]. In fact, the majority of functionals are not corrected for the self-interaction error in the exchange part as well as in the correlation part, which implies an exaggerated delocalization of the electron density that over-stabilized the transition state [55]. On the contrary, the B0KCIS, which is corrected for self-interaction error in the correlation part, gives a too positive barrier (3.6 kcal/mol) but reproducing the trend of the MP2 computation (2.5 kcal/mol). The comparison between the results given by B0LYP and B0KCIS clearly shows the very strong influence of the correlation contribution for these applications.

Proton transfer mechanisms are often of first importance in chemical or biological processes.

In spite of the apparent simplicity of this phenomenon, the theoretical description of proton transfer is a difficult problem for both post-Hartree-Fock methods and DFT approaches. In particular, it is well-known that current functionals provide very low activation energies [56][57][58]. The situation is even more involved when the motion of the proton is coupled with a rearrangement of the  electronic system, as happens in malonaldehyde. In fact, the enol form of this molecule is stabilized by an intramolecular hydrogen bond, and it is possible to switch form a stable form to the other one by a proton transfer mechanism via a transition state of C 2v symmetry [57]. (12) This spin density on the nucleus can have two origins. The first direct origin is simply provided by a unpaired electron in an orbital which has a non-zero density at the position of the nucleus. The second indirect origin is due to spin polarization of paired electrons, which have a non-zero presence probability on the nucleus, caused by the interaction with the unpaired electron under the effect of electronic correlation which acts differently depending whether the electrons have the same spin or opposite ones [START_REF] Adamo | [END_REF].

                  , * , , ) ( ) ( ) ( n n n s P P r r r
In order to test B0KCIS, we have performed computations on the methyl and vinyl radicals.

These radicals can be considered as prototype of  and  organic radicals, respectively. The computed values of hyperfine coupling constants are collected in Table VII. First of all, it is striking to notice that for this property, DFT approaches are globally more accurate than very refined correlation ab initio methods like CCSD (Coupled-Cluster Single Double) or CAS (Complete Active Space) [START_REF] Barone | V In Recent Advances in Density Functional Methods Part 1[END_REF]. For the methyl radical, the unpaired electron being in a p orbital of zero-density on the carbon atom, the hyperfine coupling constants originate exclusively from spin polarization of electrons of the 2s orbital for C or the 1s orbital for H.

The same phenomenon appears in the allyl radical. Concerning electrons of same spin, this interaction can only be treated properly with an accurate exchange functional like in PBE0 [START_REF] Adamo | [END_REF]. On the contrary, when a less accurate exchange functional is used like B, results strongly depend on error compensation with the correlation part. As a consequence, B0LYP overestimates the hyperfine coupling constant of C in the methyl radical whereas B0KCIS underestimates it.

The case of non-saturated  radicals like the vinyl radical is more complicated to treat by DFT approaches as well as by post-HF methods since it involves both direct and indirect contributions to hyperfine coupling constants [START_REF] Adamo | [END_REF]. Nevertheless, it is striking to notice the excellent value of hyperfine coupling constant of C  , dominated by the direct contribution of the unpaired electron in a  orbital of non-zero density on the nucleus, given by B0KCIS. The comparison of the results given by B0LYP and B0KCIS functionals suggests the importance of the correlation contribution for this application.

Conclusions

In this work, we have tested in a systematic way the new meta-GGA parameter-free correlation functional KCIS. Concerning atomic correlation energies, this functional is globally more accurate than other theoretical or semi-empirical correlation functionals like the correlation part of PBE or LYP, respectively. The search for the best exchange functional which can been coupled to KCIS has leaded to a new hybrid model with the Becke 88 (B) exchange functional: B0KCIS. This model, including 25 % of exact exchange and containing only one empirical parameter in the exchange part, has been tested over a varied set of physico-chemical properties and have turned out to have performances better or at least equivalent to those provided by semi-empirical exchange-correlation functionals like B3LYP.

In details, the use of B0KCIS instead of B0LYP brings a significant improvement for properties where the correlation contribution plays an important role like for atomization energies, energetic reaction barriers and magnetic properties.

While, in the present study we have searched for the best suitable exchange functional, among those already available "on the market", to be used with KCIS, we believe that it would be interesting to developed a new exchange functional entirely adapted to KCIS in order to obtain an even more accurate model. Work is in progress in such a direction.

Appendix: Expression of the KCIS correlation functional

This appendix contains the complete expression of the KCIS correlation functional. The following seven variables are involved in this meta-GGA functional:

  ,   , 2 2         g , 2 2         g ,           ,   and   :   r d F E KCIS c ) , , , , , , ( ,
                    (A.1)
where F presents the following Self-Interaction Corrected form:

) 1 , , 0 , ( ) , ( ) , , , ( ,                             g g g F GGAGAP c W GGAGAP c (A.2) with        , 2 / 1 ) 2 (             g ,         and       8 2 g W . GGAGAP c
 is the correlation energy per particle of a uniform electron gas with a gap and including some gradient corrections. It can be obtained for any spin polarization  by the formula:

                             2 , 3 , 2 , 1 , , 2 , 3 , 2 , 1 , , ) , ( ) ( ) , ( ) ( 1

Finally, the lastP

  difficult playground is represented by one of the physical constant derived from EPR spectroscopic: the isotropic hyperfine coupling constant. The isotropic hyperfine coupling constant of a nucleus n, linked to the interaction of the nuclear spin with an electronic spin by the contact Fermi mechanism, depends on the electronic spin density  are the electronic and nuclear magnetons, e g , n g the electronic and nuclear gfactors, and h is the Planck's constant. The spin density can be easily obtained from density matrices of electrons  and  , , and from the basis atomic orbitals     :

Figure 1 .Figure 2

 12 Figure 1. Mean absolute errots (MAEs) for the bond lengths of 32 molecules belonging to the G2 set. computed with the 6-311G(d.p) basis set

  energies of the unpolarized and polarized electron gas:

 10 )

 10 have been determined by a fitting procedure, starting from a model of uniform electron gas with a gap G: is the correlation energy per particle of the unpolarized uniform electron gas at 0  G , computed using the local part of the Perdew-Wang 91 (PW91) correlation functional:The gradient correction in the correlation energies per particle of the unpolarized and polarized electron gas, as to satisfy the slowly-varying limit,

  Table VI contains the most significant geometrical parameters as well as energetic barriers for the minimum and the saddle point of malonaldehyde. Whereas B0LYP over-estimates the energetic barrier with respect to the MP2 computations, (3.5 and 2.8 kcal/mol, respectively), B0KCIS gives geometrical and energetic parameters in good agreement with the reference values. Once again, the difference of 1 kcal/mol in the energetic barrier given by B0LYP and B0KCIS shows the importance of the correlation contribution and the better quality of KCIS over LYP for this kind of applications.

Table II .

 II Mean absolute errors (MAE) and maximal errors (kcal/mol) for the atomization energies of molecules belonging to the G2-1 (55 molecules) and G2 (148 molecules),

	computed with the 6-311+G(3df,2p) basis set and MP2/6-31G* geometries.
		G2-1		G2
		MAE	max.	MAE	max.
	PBE	8.2	29.1 (CO 2 )	17.2	50.5 (C 2 F 2 )
	PKZB	3.6	11.0 (O 2 )	4.5	-37.7 (SiF 4 )
	PBE0	3.1	-10.7 (SiO)	5.0	-21.7 (SiF 4 )
	B3LYP	2.4	-9.9 (SO 2 )	3.1	-20.2 (SiF 4 )
	B0LYP	5.1	-20.6 (SO 2 )	11.2	-30.1 (SiF 4 )
	PBEKCIS	11.9	35.2 (CO 2 )	24.9	66.5 (C 5 H 5 N)
	PKZBKCIS	4.1	-14.4 (C 2 H 6 )	9.7	38.8 (SiF 4 )
	PKZB opt KCIS	3.9	13.7 (O 2 )	10.0	-23.2 (C 4 H 10 )
	BKCIS	8.2	25.9 (CO 2 )	14.6	39.6 (NO 2 )
	B0KCIS	3.0	-10.5 (SO 2 )	5.3	-28.6 (SiF 4 )
	B1KCIS	3.0	-9.6 (SiO)	5.4	-28.1 (SiF4)

Table III .

 III Exact total atomic energies and differential SCF energies (Hartree), computed with

	the 6-311+G(3df,3pd) basis set.
	Atoms	Exact a	PBE	PBE0 B3LYP B0LYP PKZB	B0KCIS
	H	-0.500	0.000 -0.001	-0.002	0.002	0.004	0.002
	He	-2.904	0.014	0.011	-0.009	0.000	0.002	0.003
	Li	-7.478	0.018	0.012	-0.012	-0.003	0.004	0.002
	Be	-14.667	0.041	0.032	-0.004	0.006	0.019	0.015
	B	-24.654	0.045	0.037	-0.009	0.003	0.036	0.017
	C	-37.845	0.052	0.042	-0.012	0.000	0.056	0.015
	N	-54.589	0.060	0.049	-0.012	0.002	0.083	0.015
	O	-75.067	0.062	0.054	-0.024	-0.008	0.110	0.011
	F	-99.734	0.073	0.067	-0.028	-0.010	0.147	0.011
	Ne	-128.938 0.092	0.087	-0.022	-0.001	0.198	0.018
	MAE (H-Ne)	--	0.046	0.039	0.013	0.004	0.066	0.011
	Na	-162.255 0.099	0.085	-0.032	-0.013	0.230	0.006
	Mg	-200.053 0.108	0.092	-0.040	-0.023	0.267	0.001
	Al	-242.346 0.119	0.098	-0.041	-0.024	0.308	-0.001
	Si	-289.359 0.134	0.108	-0.035	-0.020	0.351	-0.001
	P	-341.259 0.155	0.124	-0.023	-0.008	0.401	0.005
	S	-398.110 0.170	0.135	-0.024	-0.010	0.446	0.005
	Cl	-460.148 0.189	0.150	-0.020	-0.006	0.495	0.007
	Ar	-527.540 0.209	0.165	-0.014	0.001	0.545	0.007
	MAE (Na-Ar)	--	0.148	0.120	0.029	0.013	0.380	0.004
	Total MAE	--	0.091	0.075	0.020	0.008	0.206	0.008
	a) ref. 35

Table IV .

 IV Dissociation energies (kcal/mol) and main geometrical parameters (Å) for H-bond dimers. The energies (corrected for BSSE error) have been computed with the 6-

	311++G(3df,3pd) basis set with the 6-311G(d,p) geometries
		B3LYP	B0LYP	PKZB	B0KCIS	Best	Exp. b
						ab initio a
	(H 2 O) 2
	d(OO)	2.888	2.892	2.883	2.912	2.925	2.952
	d(H .. O)	1.926	1.931	2.364	1.952
	E	4.6	4.5	1.9	3.8	4.7	5.4  0.7
	(HF) 2
	d(FF)	2.707	2.710	2.888	2.734	2.76	2.73/2.74
	d(H .. F)	1.780	1.786	1.956	1.818
	E	4.7	4.7	3.0	3.8	4.6	4.2
	(HCl) 2
	d(ClCl)	3.917	3.947	4.387	4.009	3.790	3.75/3.79
	d(H .. Cl)	2.626	2.657	3.093	2.724
	E	1.2	0.9	0.8	0.8	2.0	1.4/2.0
	a) ref. 43, 45, 46; b) ref. 44, 47, 48.

Table V .

 V Complexation energies of the ion-molecule complex (E comp ), activation energies (E # ) and energetic barriers relative to reagents (E r ) for the considered S N 2 reaction, computed at different levels of approximation. All the energies (kcal/mol) have been

	computed with the 6-311++G(3df,3pd) base set and the 6-311+G(d,p) geometries.
		G2+(MP2)	B3LYP B0LYP PKZB	B0KCIS	Exp.
	 E comp E # E r	10.6 13.0 2.5	9.8 9.1 -0.7	9.7 9.7 0.0	9.6 9.6 0.0	8.4 11.9 3.6	12.2  2 13.3  2 1  1

Table VI .

 VI Main geometrical parameters (Å) of the H-bond and relative energies (kcal/mol) for the minimum (M) and the saddle point (SP) of malonaldehyde. The energies have been computed with the 6-311++G(3df.3pd) basis set with the 6-311G(d.p) geometries.

		MP2		B3LYP	B0LYP	B0KCIS
		M	SP	M	SP	M	SP	M	SP
	d(OO)	2.581 2.355 2.580 2.367 2.593 2.363	2.564	2.356
	d(H .. O)	1.678 1.197 1.687 1.208 1.711 1.207	1.666	1.202
	d(OH)	0.991 1.197 0.999 1.208 0.992 1.207	0.996	1.202
	E	0.0	2.8	0.0	3.0	0.0	3.5	0.0	2.5

Table VII .

 VII Isotropic hyperfine coupling constants (a. Gauss) for the methyl and vinyl radical computed with the EPR-III basis set.

		B3LYP B0LYP B0KCIS CCSD(T)	Exp.
	methyl a(C) 29.8	31.7	24.9	27.8	28.4
		a(H) -23.3 -24.3 -26.2	-24.6	-25.1
	vinyl a(C a ) 113.1 112.7 107.6	121	107.6
		a(C b ) -5.8	-5.8	-5.2	-8	-8.6
		a(H a ) 17.8	17.7	12.3	9	13.3
		a(H anti ) 41.9	41.9	40.7	28	34.2
		a(H syn ) 65.6	65.5	62.2	48	68.5
	MAE	3.8	3.9	3.3	6.6	--