Diego Cattaruzza
email: cattaruzza@emse.fr

Nabil Absi
email: absi@emse.fr

Dominique Feillet
email: feillet@emse.fr

Olivier Guyon
email: guyon@emse.fr

Xavier Libeaut
email: xavier.libeaut@uco.fr

The Multi Trip Vehicle Routing Problem with Time Windows and Release Dates

  L'archive ouverte pluridisciplinaire

Introduction

The well known Vehicle Routing Problem (VRP) is an NP-hard combinatorial optimization problem where a set of geographically scattered customers has to be served by a fleet of vehicles minimizing routing costs and respecting capacity constraints on vehicles.

In some city distribution systems, goods are delivered to city distribution centers (CDC) and then delivered to final customers by vehicles with limited capacity. They can be back to the depot much earlier than the end of the working day and then be available for another delivery trip. This introduces the multi trip aspect.

Customers usually ask to be served within a certain time interval. Meeting these intervals is vital for the carrier: delays mean losing reliability and trustworthiness and often penalties need to be payed. Then, time windows should be considered and associated with each customer.

Finally, merchandise can be delivered to the CDC all day long. This means that they are not completely available at the CDC at the beginning of the planning horizon. The concept of release date is then associated with each merchandise.

In this paper we address the Multi Trip Vehicle Routing Problem with Time Windows and Release Dates (MTVRPTWR). It is noteworthy that the problem is deterministic even if the merchandise continuously arrives to the CDC during the day: the release dates are supposed to be known before the working day starts.

The paper is structured as follows. Section 2 reviews the literature, Section 3 formally defines the problem. In Section 4 the hybrid genetic algorithm (HGA) we proposed to solve the MTVRPTWR is explained. Results are reported in Section 5, while Section 6 concludes the paper.

Literature Review

In multi-echelon distribution systems, delivery to customer is not direct, but goods are firstly delivered to CDC. Final distribution trips need to take into account the moment goods become available at the CDC. Normally, only the availability of the vehicles to perform deliveries is required. When goods must be loaded immediately after unloading, synchronization aspects can be introduced (Drexl [6]).

We model time-interdependency among vehicles introducing the concept of release date associated with merchandise. To the best of our knowledge no previous work is done on it.

Few papers study the Multi Trip VRP with Time Windows (MTVRPTW) and most of all use exact methods to solve it.

Azi et al. [START_REF] Azi | An exact algorithm for a single-vehicle routing problem with time windows and multiple routes[END_REF] propose an exact algorithm to solve the single vehicle MTVRPTW. The solution approach exploits an elementary shortest path algorithm with resource constraints. In the first phase all non-dominated paths are calculated. Then the shortest path algorithm is applied to a modified graph. Each node represents a non-dominated route and two nodes are connected by an arc whether it is possible to serve the two routes consequently and they do not serve common customers. Solomon instances are used with different values of time horizon. 16 instances out of 54 with 100 customers are solved to optimality. Azi et al. [START_REF] Azi | An exact algorithm for a vehicle routing problem with time windows and multiple routes[END_REF] address the MTVRPTW with a column generation approach embedded within a branch-and-price algorithm. A set packing formulation is given for the master problem and each column represents a working day. Since each pricing problem is an elementary shortest path with resource constraints, a similar approach to the one proposed in Azi et al. [START_REF] Azi | An exact algorithm for a single-vehicle routing problem with time windows and multiple routes[END_REF] is applied. As in Azi et al. [START_REF] Azi | An exact algorithm for a single-vehicle routing problem with time windows and multiple routes[END_REF], Solomon instances are considered and a time horizon is introduced. Due to the limitation of the algorithm, the authors focus on instances formed by the first 25 or 40 customers of each Solomon's instance. The algorithm can also solve few instances of size 50. 24 instances out of 27 of size 25 are solved within 30 hours. When the number of customers increases to 40, only 9 instances out of 27 can be solved to optimality. Hernandez et al. [START_REF] Hernandez | A new exact algorithm to solve the multitrip vehicle routing problem with time windows and limited duration[END_REF] use a similar approach of Azi et al. [START_REF] Azi | An exact algorithm for a vehicle routing problem with time windows and multiple routes[END_REF]. A set covering formulation is given for the problem and each column represents a trip instead of a working day. Results on the same instances of Azi et al. [START_REF] Azi | An exact algorithm for a vehicle routing problem with time windows and multiple routes[END_REF] show that 25 out of 27 instances with 25 customers and large time horizon and 2 vehicles and 22 out of 27 instances of size 50 with large time horizon and 4 vehicles have been solved within 30 hours. All the three exact methods cited use the algorithm proposed by Feillet et al. [START_REF] Feillet | An exact algorithm for the elementary shortest path problem with resource constraints: application to some vehicle routing problems[END_REF] for the elementary shortest path problem.

Battarra et al. [START_REF] Battarra | An adaptive guidance approach for the heuristic solution of a minimum multiple trip vehicle routing problem[END_REF] study an extension of the MTVRPTW where products are clustered in different commodities that cannot be transported in the same vehicle. The approach they propose is a two-stage sequential heuristic: they generate a set of feasible routes considering each commodity independently, then, routes are assigned to vehicles.

Problem definition and notation

The MTVRPTWR can be defined on a completed undirected graph G = (V, E), where V = {0, . . . , N } is the set of vertices and E = {(i, j)|i, j ∈ V, i < j} the set of edges. Vertex 0 represents the depot, where a fleet of M identical vehicles with capacity Q is based. Vertices 1, . . . , N represent the customers. With each customer is associated a demand D i and a Time Window (TW) [E i , L i ] during which the service should start and a service time S i . Arriving at customer location before E i is allowed. Since the service cannot start earlier than E i the driver must wait. Late arrivals at customer location are forbidden. Moreover, the quantity D i of product requested by customer i is available at the depot not earlier than R i . R i is called release date. For brevity, we say R i is associated with customer i, instead of R i is associated with the demand D i requested by customer i. It is possible to travel from i to j incurring in a travel time T ij . A vehicle, before performing each trip, needs to wait until goods it is going to deliver are available at the depot, i.e., it cannot start the trip before the maximal release date associated with customers it has to serve.

A time horizon T H is given and establishes the duration of the working day. It must be respected and it can be viewed as a TW associated with the depot. Thus, it is assumed that

[E 0 , L 0 ] = [0, T H ]. Moreover D 0 = 0.
The MTVRPTWR calls for the determination of a set of routes and an assignment of each route to a vehicle, such that the total travel time is minimized and the following conditions are satisfied: (a) each route starts and ends at the depot;

(b) each route starts not earlier than the largest R i associated with customers assigned to the route itself; (e) the sum of the demands of the customers in any route does not exceed Q;

(f) each vehicle arrives at the depot after performing its last route not later then T H (TW associated with the depot is respected).

It is also supposed that each customer i could be served by a return trip, i.e, R i + T 0i + T i0 ≤ T H and D i ≤ Q (otherwise no feasible solution would exist).

The following notation is used in the remaining. The symbol σ indicates a route. The capital Σ indicates the set of routes assigned to a vehicle. For easiness, it is called journey. Then, a journey is formed by different routes. The symbol ⊕ indicates the concatenation of paths (partial routes) or routes. For example (v 1 , . . . , v n ) ⊕ (w 1 , . . . , w m ) is the concatenation of two paths (that results in a route if v 1 and w m are the depot). σ 1 ⊕ σ 2 means that route σ 1 is performed right before σ 2 on the same vehicle. The time a vehicle needs to visit customers on route σ i is indicated by τ σ i , while the time the vehicle leaves the depot to perform σ i is indicated by T σ i .

The symbol "∈" describes "belonging". For example σ ∈ Σ means route σ belongs to journey Σ, v ∈ σ means customer v belongs to (i.e., it is served by) route σ, and so on.

It is noteworthy that in absence of release dates it would hold T σ i+1 = T σ i + τ σ i . On the other hand, considering release dates it holds

T σ i+1 = max {T σ i + τ σ i , max v i ∈σ i {R v i }}.
A customer v i is feasible if the vehicle arrives at its location before L i , infeasible otherwise. A route is feasible if it serves only feasible customers, infeasible otherwise. A journey is feasible if it is composed by feasible routes, infeasible otherwise.

Method

We propose a genetic algorithm (HGA) for the MTVRPTWR. The main scheme of our approach is illustrated in Algorithm 1. The algorithm is an adaptation of the procedure proposed by Cattaruzza et al. [START_REF] Cattaruzza | A memetic algorithm for the multi trip vehicle routing problem[END_REF] Algorithm 1 Hybrid Gentic Algorithm sketch 1: Initialize population 2: while Termination criteria are not met do 3:

Select parent chromosomes Ψ P 1 and Ψ P 2

4:

Generate a child Ψ C

5:

Improve

Ψ C with LS 6: if Ψ C is infeasible then 7: Repair Ψ C 8:
end if

9:
Insert Ψ C in the population 10:

if Dimension of the population is bigger or equal than π + µ then 11:

Select survivors 12:

end if 13: end while (we refer to this paper for more details). First, π random chromosomes are generated in order to create an initial population. These individuals are then improved with Local Search (LS). The classical binary tournament procedure is used for selection and children are created by crossing the selected parents (with the OX operator). Infeasible chromosomes are repaired (Section 4.4). When the population reached the dimension of π + µ, µ individuals are eliminated accordingly to their quality and their diversification contribution to the population (as proposed by Vidal et al. [START_REF] Vidal | A hybrid genetic algorithm for multi-depot and periodic vehicle routing problems[END_REF]).

In the remaining, Sections 4.2 and 4.3 provide details for the LS while Section 4.5 describes AdSplit procedure used to obtain MTVRPTWR solutions from a chromosome.

Solution representation and search space

A chromosome is a sequence (permutation) Ψ = (Ψ 1 , . . . , Ψ N ) of the N client nodes, without trip delimiters. Ψ can be viewed as a TSP solution that has to be turned in a feasible MTVRPTWR solution by splitting the chromosome (inserting trip delimiters and assigning trips to vehicles). From that point of view, Ψ is usually called a giant tour. From a giant tour Ψ, different MTVRPTWR solutions can be constructed depending on the way Ψ is split.

During the search phase, overload and TW violations are allowed and penalized in the fitness function respectively by coefficients λ and θ.

The procedure AdSplit (see Section 4.5) is used to get a MTVRPTWR solution ξ from Ψ. The two following notations are introduced: T Σ (ξ) and T W Σ (ξ) are respectively the travel time and the TW violation of journey Σ in solution ξ. L σ (ξ) is the load of trip σ where σ ∈ Σ is a route σ of journey Σ. The fitness F (Ψ) of the chromosome Ψ is the cost of the best solution ξ found by AdSplit and it is defined as follows:

F (Ψ) = c(ξ) = M Σ=1 T Σ (ξ) + θ M Σ=1 T W Σ (ξ) + λ M Σ=1 σ∈Σ max{0, L σ (ξ) -Q} (1)
When confusion cannot arise, solution ξ will be omitted in the notation. The chromosome Ψ is called feasible (infeasible) if AdSplit obtains, from Ψ, a feasible (infeasible) solution ξ.

Local Search for VRP with Time Windows

LS in presence of TW becomes more complicated than in the classic VRP. In particular, feasibility cannot be checked in constant time straightforwardly. In our approach, we adapt the scheme proposed by Vidal et al. [START_REF] Vidal | A hybrid genetic algorithm with adaptive diversity management for a large class of vehicle routing problems with time windows[END_REF] (that is in turn an extension of the scheme proposed by Nagata et al. [START_REF] Nagata | A penalty-based edge assembly memetic algorithm for the vehicle routing problem with time windows[END_REF]) for a large class of problems with TWs. Given a path ρ, the quantities T (ρ), T W (ρ), E(ρ), L(ρ), D(ρ) respectively denote the minimum duration, the minimum TW violation of ρ, the earliest and the latest date service can start at the first customer of ρ (that can be the depot) allowing minimum duration and TW penalization, and the cumulative demand of served customer. Given two paths ρ 1 and ρ 2 , the following relations hold:

T (ρ 1 ⊕ ρ 2 ) = T (ρ 1 ) + T (ρ 2 ) + T v 1 ,v 2 + ∆ W T ; (2) T W (ρ 1 ⊕ ρ 2 ) = T W (ρ 1 ) + T W (ρ 2 ) + ∆ T W ; (3) E(ρ 1 ⊕ ρ 2 ) = max {E(ρ 2 ) -∆, E(ρ 1 )} -∆ W T ; (4) L(ρ 1 ⊕ ρ 2 ) = min {L(ρ 2 ) -∆, L(ρ 1 )} + ∆ T W ; (5) D(ρ 1 ⊕ ρ 2 ) = D(ρ 1 ) + D(ρ 2 ); (6)
where

∆ = T (ρ 1 ) -T W (ρ 1 ) + T v 1 ,v 2 ; ∆ W T = max {E(ρ 2 ) -∆ -L(ρ 1 ), 0} ; ∆ T W = max {E(ρ 1 ) + ∆ -L(ρ 2 ), 0} ; v 1 last customer in ρ 1 and v 2 first customer in ρ 2 .
Knowing those quantities allows to evaluate classical LS operators in constant time. The quantities defined above for all the paths and their reverse (needed, for example, to evaluate 2-opt moves) can be calculated in a preprocessing phase.

Local search: introduction of release dates and application to the multi-trip case

The problem we address has two more characteristics than the VRPTW: vehicles can perform several trips and goods become available for final distribution throughout the time horizon. The consequence is Singapore, August 4-8, 2013 MIC 2013: The X Metaheuristics International Conference id-5 that vehicles can start a trip at t > 0, because the trip they are going to accomplish is not the first, or either because they wait for goods to be available.

Vidal et al. [START_REF] Vidal | A hybrid genetic algorithm with adaptive diversity management for a large class of vehicle routing problems with time windows[END_REF] prove relations 2-6 by induction on the concatenation operator proving the following stronger relations:

T (ρ)(t) = T (ρ) + max {0, E(ρ) -t} ; (7) T W (ρ)(t) = T W (ρ) + max {0, t -L(ρ)} . ( 8
)
Relation 8 is used to calculate effective TW violation due to later departing from the depot. We introduce R(ρ) as the largest release date of customers served in ρ. We define

R(v i ) = R v i and it holds R(ρ 1 ⊕ ρ 2 ) = max {R(ρ 1 ), R(ρ 2 )} (9)
Thus, release dates, do not introduce any complexity into the moves evaluation scheme introduced in Section 4.2.

It remains to consider the multi-trip aspect. The idea is to consider each route σ i as it would start at T σ i = 0 and calculate variation in TW violations and time completion based on the effective starting time of the trip itself, i.e., considering 8can be used to calculate T W (σ i ) for all the trips in Σ, while T σ i+1 can be calculated as

T σ i = max R(σ i ), T σ i-1 + τ σ i-1 . Given a journey Σ = (σ 1 ⊕ • • • ⊕ σ n ), Equation
t =      E(σ i ) + T (σ i ) -T W (σ i ) if T σ i ≤ E(σ i ), T σ i + T (σ i ) -T W (σ i ) if E(σ i ) < T σ i ≤ L(σ i ), L(σ i ) + T (σ i ) -T W (σ i ) if L(σ i ) < T σ i . (10)
T σ i+1 = max{t, R(σ i+1 )} (11)
The effect of a move on trip σ i ∈ Σ can be propagated to each of the following routes by means of a variation in completion time and in TW penalization. Using Equations 7 and 8 this variation can be locally calculated in constant time. Complete calculation requires O(Σ max ), where Σ max indicates the maximum numbers of trips among all journeys.

Repair Procedure

Infeasible chromosomes are repaired applying LS with penalization factors multiplied by 10 depending on the nature of the infeasibility. If the resulting chromosome is still infeasible, λ and/or θ are multiplied again by 10 and LS is re-applied.

A split algorithm for the MTVRPTWR 4.5.1 Auxiliary graph construction

The splitting procedure AdSplit, is an adaptation of the procedure proposed by Prins in [START_REF] Prins | A simple and effective evolutionary algorithm for the vehicle routing problem[END_REF]. It works on an auxiliary graph H = (V ′ , A ′ ). V ′ contains N + 1 nodes indexed from 0 to N . Arc (i, j), i < j, represents a trip σ j i+1 serving customers from Ψ i+1 to Ψ j in the order they are in a given chromosome Ψ, i.e., σ j i+1 = (0, Ψ i+1 , . . . , Ψ j , 0). Since we are looking for feasible solutions, we could consider only feasible trips with the corresponding routing cost c ij associated with arc (i, j), i.e.,

c ij = T (σ j i+1 )
The best solution associated with chromosome Ψ would be represented by the path that goes from node zero to the node Ψ N respecting time and capacity constraints.

In fact, we do not have the guarantee that such a feasible path (and then such a feasible solution) exists. We, then, allow infeasibility with respect to time and load constraints and we associate with each arc (i, j) the following cost

c ij = T (σ j i+1 ) + θT W (σ j i+1 ) + λ max{0, D(σ j i+1 ) -Q}. ( 12
)
It is noteworthy that the arc costs do not take into account the position of the trip in the journey, but it is the (penalized) cost of the trip as it would be performed as first (i.e., at t = 0). Therefore, the (contingent) TW violation due to later departure is not taken into account. As a consequence, considering the arcs that form the shortest path on H to construct the solution (as in the original procedure in Prins [START_REF] Prins | A simple and effective evolutionary algorithm for the vehicle routing problem[END_REF] for the VRP and in Cattaruzza et al. [START_REF] Cattaruzza | A memetic algorithm for the multi trip vehicle routing problem[END_REF] for the MTVRP) could provide a solution with a cost much higher than the best one. We then propose a labelling procedure (explained in Section 4.5.2) that obtains a solution for the MTVRPTWR starting from the graph H constructed on Ψ.

Assignment of trips to vehicles

In the MTVRPTWR context in particular and in the MTVRP context in general, more than one trip can be assigned to the same vehicle. TW penalization deeply depends on the time routes leave the depot and this aspect cannot be considered by the constant costs associated with arcs on H (Section 4.5.1).

We propose the following labelling procedure that enumerates all the possible paths from node 0 to N in Ψ and construct the best solution associated with the label with lower cost on node N .

Starting from node 0, labels are progressively extended along the graph. Each label L has M + 3 fields: the first M fields store vehicle travel times in decreasing order, the (M + 1) th field memorizes the total load infeasibility, the (M + 2) th the predecessor node, and the last field keeps the cost of the partial solution evaluated using Equation 1and equivalent to the cost c(L) of label L. When extending a label, M new labels are constructed, one for each possible allocation of the new trip to a vehicle.

When node N is reached, the label L with minimum cost c(L) associated with node N is selected and the related solution is constructed (going backwardly trough the graph).

To speed up the procedure, dominated labels, accordingly with the following dominance rule, are discarded: let L 1 and L 2 be two labels associated with the same node i ∈ V ′ . L 1 strongly dominates L 2 if and only if

c(L 1 ) + θ M j=1 δ j (L 1 , L 2 ) ≤ c(L 2 ) (13)
where c(L) is the cost associated with label L, δ j (L 1 , L 2 ) = max{0, T j (L 1 ) -T j (L 2 )} and T j (L) is the (partial) travel time of vehicle j associated with label L. Roughly speaking, given two labels L 1 and L 2 , extending L 1 TW violation is penalized as much as possible while it is not extending L 2 in the same way. If Inequality 13 holds, L 2 cannot be extended in a better way than L 1 , then, L 2 is eliminated. Preliminary computational experiments (Section 5) show the time inefficiency of the procedure. For this reason a heuristic version of the dominance rule is proposed and a parameter γ ≥ 1 is introduced. L 1 weakly dominates L 2 if and only if

c(L 1 ) + θ M j=1 δ j (L 1 , L 2 ) ≤ γc(L 2 ), (14)
c(L 1 ) ≤ c(L 2 ). (15)
It is noteworthy that for γ = 1, the weak dominance rule is equivalent to the strong version. When γ > 1, Inequality 14 is easily satisfied and a larger number of labels can be eliminated. Condition 15 is added because when γ > 1 label L 2 can be dominated by label L 1 even if c(L 2 ) < c(L 1 ). This cannot happen if γ = 1. Using the weak relation one expects the solution to be obtained quicker. On the other side the best decomposition of the chromosome can be missed.

Singapore, August 4-8, 2013

MIC 2013: The X Metaheuristics International Conference id-7

The value of γ is dynamically adapted during the process accordingly with the number of the label associated with each node. Precisely the following scheme is adopted:

γ = γ + |L i | 1000L threshold if |L i | > L threshold γ -L threshold 1000|L i | if |L i | < L threshold (16)
where |L i | is the number of labels associated with node i and L threshold is a threshold parameter that indicates the number of labels that should be kept associated with each node. If after Equation 16 is applied, γ is lower than 1, it is set back to 1. This is natural, since for γ < 1 dominated labels would be kept. The lower L threshold the quicker is the procedure and the poorer is the quality of the solution obtained.

Discussion

In this section the computational results obtained are discussed. First of all we explain how we generated a set of 19 instances for the MTVRPTWR, since to the best of our knowledge there is not a benchmark of instances available in the literature for the problem. The value of the L threshold is determined (Section 5.2). Finally results are presented (Sections 5.3-5.4).

Instances generation

We define a generator of instances for the MTVRPTWR. The time horizon T H is randomly drawn from [500, 520] and is represented by a TW [0, T H ] associated with the depot. The capacity Q of each vehicle is fixed to 100. The service time at the depot S 0 is 20, while each customer i requires a service time S i = 5. Clients and the depot are located on a [0, 50] × [0, 50] squared Cartesian system portion. The depot is centrally located. For each customer, its abscissa and ordinate are randomly drawn in the interval [0, 50]. Given two customers i and j, the travel time T ij is the rounded Euclidean distance between them. Customer's request is randomly drawn in the interval [START_REF] Clarke | Scheduling of vehicles from a central depot to a number of delivery points[END_REF]15]. The Clarke and Wright heuristic (Clarke and Wright [START_REF] Clarke | Scheduling of vehicles from a central depot to a number of delivery points[END_REF]) is then applied to these data and a VRP solution is obtained. Trips σ obtained are concatenated as long as the journey respects depot's TW. Corresponding T σ are calculated. If a trip cannot be assigned to any vehicle, a new journey is initialized with it. M is fixed to the number of vehicles which are needed to assign each trip to a journey. For each customer i, two dates are calculated: the starting of the service s i time when the vehicle leaves the depot at t = 0 and s i when the vehicle leaves the depot as late as possible respecting T H . Then s i is calculated as the average of s i and s i . L i is then calculated as follows:

• If s i < 240

-L i = 240 with a probability of 0.7; -L i = 360 with a probability of 0.2; -L i = 480 with a probability of 0.1.

• If 240 ≤ s i < 360
-L i = 360 with a probability of 0.6; -L i = 480 with a probability of 0.4.

• Otherwise L i = 480.

The release dates R i are determined as follows:

• R i = 0 with a probability of 0.5.

• Otherwise R i is randomly drawn from [0, R σ ], where R σ = min j∈σ {T σs j , L js j } and i ∈ σ.

Finally E i = R i + S 0 + T 0i . A

Determination of L threshold

To determine the value of L threshold , 60 chromosomes for an instance with N = 150 and M = 6 are randomly generated. Then, they are evaluated by the AdSplit procedure with the use of the strong dominance rule and the weak dominance rule with different values of L threshold . Results are reported in Table 1. It can be noticed that using the strong dominance rule more than 3 minutes are averagely needed to split a chromosome. This justifies the introduction of the weak dominance rule. Still from Table 1 it can be noticed that decreasing L threshold , the procedure is quicker but the quality of the obtained solutions is slightly reduced. It has been decided to set L threshold = 5.

Comparison with a classical Tabu Search algorithm

The procedure is run once over all instances for 5 minutes. Results are compared with those obtained by a Tabu Search (TS) algorithm that we developed to obtain reference values. The TS makes use of the classical neighborhoods insert and swap. Seven structures are considered combining them differently.

The TS is run 6 times on each instance with each one of the different structures for a total of 42 times (each stopped after 5 minutes of computation). Results are reported in

Comparison with Hernandez et al. [9]

Hernandez et al. [START_REF] Hernandez | An exact algorithm to solve the multi-trip vehicle routing problem with time windows[END_REF] propose an exact method for the MTVRPTW and are able to solve 30 out of 54 instances they construct from Solomon's instances for the VRPTW. We run our algorithm on those instances (R i = 0 for all i). Results are reported in Table 3.

The first column reports the instance name ON M N , where ON is Solomon's instance original name. The second column reports the optimal value, the third the cost of the solution found by HGA, while the last column reports the CPU time (in seconds). The procedure is stopped when the optimal solution is found or after 5 minutes. HGA finds the optimal value 25 out of 30 times missing it just once for instances with M = 2 and N = 25. The average gap from the optimal value is 0.4%.

  (c) each customer is visited by exactly one route; Singapore, August 4-8, 2013 MIC 2013: The X Metaheuristics International Conference id-3 (d) service at customer i starts during the associated TW [E i , L i ];

Table 1 :

 1 set of 19 instances is then created with 20 to 200 customers and 1 to 7 vehicles. Comparison of dominance rules: average results obtained by AdSplit procedure on 60 chromosome for an instance with M = 6 and N = 150 (times in seconds)

	Strong						Weak
			L threshold = 20	L threshold = 10	L threshold = 5	L threshold = 1
	Cost	Time	Cost		Time	Cost	Time	Cost	Time	Cost	Time
	1138.94 193.6 0.13% -99.74% 0.19% -99.87% 0.21% -99.94% 7.67% -99.99%
					0.51		0.23		0.09		0.01
			M N Best TS Average TS HGA Best known
					cost	% gap	% gap	cost
			1	20	203	0.0		0.0	203
			2	30	309	0.0		0.0	309
			2	40	347	0.0		0.0	347
			2	50	384	2.9		0.0	384
			3	60	491	1.6		0.0	491
			3	70	469	1.9		0.0	469
			3	80	571	1.9		-0.4	569
			4	90	570	5.8		0.9	570
			4 100	642	4.5		-0.2	641
			4 110	786	2.3		-1.9	771
			5 120	755	5.2		-0.7	750
			5 130	827	3.1		-1.9	811
			5 140	952	4.0		-4.5	909
			6 150	971	7.5		-0.9	962
			6 160	1028	5.1		-1.9	1008
			6 170	1033	3.1		2.0	1033
			7 180	1197	5.1		-7.5	1107
			7 190	1208	7.4		-3.6	1164
			7 200	1245	6.7		-2.6	1213

Table 2 :

 2 Results on the instances for the MTVRPTWR

Table 2 .

 2 The first two columns report the number of vehicles and customers on each instance. Column Best TS reports the cost of the best solution obtained by the TS procedure over the 42 runs, while column Average TS reports the gap of the average cost of solution obtained over the 42 runs with respect to the Singapore, August 4-8, 2013

	Instance	opt	HGA CPU Time
	C201 2 25	380.8	387.7	300.0
	C202 2 25	368.6	368.6	2.1
	C203 2 25	361.7	361.7	0.3
	C204 2 25	358.8	358.8	0.0
	C205 2 25	377.2	377.2	0.0
	C206 2 25	367.2	367.2	0.2
	C207 2 25	359.1	359.1	0.3
	C208 2 25	360.9	360.9	15.5
	R201 2 25	554.6	554.6	0.1
	R202 2 25	485.0	485.0	1.2
	R203 2 25	444.2	444.2	0.3
	R204 2 25	407.5	407.5	8.0
	R205 2 25	448.4	448.4	0.0
	R206 2 25	413.9	413.9	0.5
	R207 2 25	400.1	400.1	0.3
	R208 2 25	394.3	394.3	2.7
	R209 2 25	418.3	418.3	0.1
	R210 2 25	448.3	448.3	117.4
	R211 2 25	400.1	400.1	0.2
	RC201 2 25 660.0	660.0	27.8
	RC202 2 25 596.8	596.8	1.2
	RC203 2 25 530.1	530.1	0.0
	RC205 2 25 605.3	605.3	0.9
	RC206 2 25 575.1	575.1	27.1
	RC207 2 25 528.2	528.2	0.1
	R201 4 50	909.8	922.7	300.0
	R202 4 50	816.0	816.0	69.6
	R205 4 50	807.3	816.2	300.0
	RC201 4 50 1096.6 1097.6	300.0
	RC202 4 50 1001.6 1038.6	300.0

Table 3 :

 3 Comparison with results of Hernandez et al.[START_REF] Hernandez | An exact algorithm to solve the multi-trip vehicle routing problem with time windows[END_REF] value in Best TS. Column HGA reports the cost of the solution obtained by the HGA procedure run just once. Column Best known reports the cost of the best known solution. HGA is better than TS on 11 instances and provides equivalent solutions on 6 instances.

Singapore, August 4-8, 2013

Conclusion

In this paper, we introduced the Multi Trip Vehicle Routing Problem With Time Windows and Release Dates. This problem is particularly interesting in the City Logistics context. We proposed a hybrid GA based on a new split procedure to face it. A set of instances has been introduced. Comparison with optimal solution of Hernandez et al. [START_REF] Hernandez | An exact algorithm to solve the multi-trip vehicle routing problem with time windows[END_REF] proves the efficiency of our algorithm.


Julien Toulouse
email: toulouse@lct.jussieu.fr

Andreas Savin
email: savin@lct.jussieu.fr

Local density approximation for long-range or for short-range energy functionals?

Density functional methods were developed, in which the Coulomb electron-electron interaction is split into a long-and a short-range part. In such methods, one term is calculated using traditional density functional approximations, like the local density approximation. The present paper tries to shed some light upon the best way to do it by comparing the accuracy of the local density approximation with accurate results for the He atom.

I. INTRODUCTION

In recent years, there has been a growing interest in approaches of density functional theory (DFT) [1] based on a long-range/short-range decomposition of the Coulomb electron-electron interaction. The idea is to use different, appropriate approximations for the long-range and the short-range contributions to the exchange and/or correlation energy density functionals of the Kohn-Sham (KS) scheme [2].

Various schemes combining a (semi)local short-range exchange energy functional approximation with an "exact" long-range exchange energy expression have been proposed (see, e.g., Refs. [3][4][5][6][7][8], which allow to impose the exact Coulombic asymptotic behavior 1/r of the exchange interaction. This turned out to be important for charge transfer, van der Waals interactions, etc. But opposite approaches combining a (semi)local long-range exchange functional approximation with an "exact" (or "hybrid") short-range exchange have also been used (see, e.g., Refs. 9-13), which allow to introduce exact exchange in solid-state calculations without the computationallydemanding exact long-range contribution. There is also a physical reason for doing it: long-range correlations are not well treated by (semi)local density approximations. Treating only the exchange exactly destroys the balance of errors, which is important for metals, gaps, etc. Contemplating this two opposite approaches, one can ask: is it preferable to use (semi)local density functional approximations for the long-range or for the short-range contribution to the exchange energy?

Nozières and Pines [14] have first used the idea of decomposing the correlation energy of the uniform electron gas into long-range and short-range contributions to facilitate its calculation. In the context of DFT, a few schemes combining a (semi)local short-range density functional approximation with a long-range correlation energy calculated by other means have been proposed (see, e.g., Refs. 15-19) to deals with for example neardegeneracy or long-range van der Waals interactions. In-deed, it is well known that (semi)local density functional approximations are appropriate for the short-range contribution (see, e.g., Refs. 20,21). However, for a given decomposition of the Coulomb interaction, there are at least two possible definitions for a short-range correlation functional. It can be defined either as the difference between the Coulombic correlation functional and a longrange correlation functional associated to the long-range part of interaction (as in Ref. 17), or directly from the short-range part of the interaction (as in Ref. 16). Therefore, one can wonder: what is the preferable definition for a short-range correlation functional as regards as the accuracy of (semi)local density functional approximations?

This work sheds some light on these two questions. Taking the example of the He atom, we test the accuracy of the local density approximation (LDA) to the long-range and short-range contributions to the exchange and correlation energies for a given decomposition of the interaction. Atomic units (a.u.) are used throughout this paper.

II. LONG-RANGE AND SHORT-RANGE DENSITY FUNCTIONALS

We define in this section the long-range and shortrange density functionals that we consider in this work.

Let's define first a general exchange functional for an arbitrary electron-electron interaction w ee (r)

E x [n; w ee ] = Φ[n]| Ŵee |Φ[n] -E H [n; w ee ], (1)
where Ŵee = (1/2) n2 (r 1 , r 2 )w ee (r 12 )dr 1 dr 2 is the interaction operator, expressed with the pair-density operator n2 (r 1 , r 2 ), Φ[n] is the KS determinant and E H [n; w ee ] = (1/2) n(r 1 )n(r 2 )w ee (r 12 )dr 1 dr 2 is the Hartree functional. Likewise, the general correlation functional associated to w ee (r) writes

E c [n; w ee ] = min Ψ→n Ψ| T + Ŵee |Ψ -Φ[n]| T + Ŵee |Φ[n] , (2)
where the constrained-search formalism [22] has been used.

We use the following decomposition of the Coulomb electron-electron interaction w coul ee (r) = 1/r w coul ee (r) = w lr,µ ee (r) + w sr,µ ee (r),

where w lr,µ ee (r) = erf(µr)/r is a long-range interaction and w sr,µ ee (r) = erfc(µr)/r is its short-range complement. The parameter µ which controls the range of the decomposition varies between 0 and ∞. For µ = 0, the long-range interaction vanishes, w lr,µ=0 ee = 0, and the short-range interaction reduces to the Coulomb interaction w sr,µ=0 ee = w coul ee . In the limit µ → ∞, the longrange interaction reduces to the Coulomb interaction w lr,µ→∞ ee = w coul ee , and the short-range interaction vanishes, w sr,µ→∞

ee = 0. The Coulombic exchange functional E coul x [n] = E x [n; w coul ee ] is consequently decomposed as E coul x [n] = E lr,µ x [n] + E sr,µ x [n], (4)
where

E lr,µ x [n] = E x [n; w lr,µ
ee ] is the long-range exchange functional associated to the interaction w lr,µ ee , and

E sr,µ x [n] = E x [n; w sr,µ
ee ] is the short-range exchange functional associated to the interaction w sr,µ ee . One can also define a long-range correlation functional associated to the interaction w lr,µ ee , E lr,µ c [n] = E c [n; w lr,µ ee ], and a short-range correlation functional associated to the interaction w sr,µ ee , E sr,µ c

[n] = E c [n; w sr,µ ee ]. However, because E c [n; w ee ] is not linear with respect to w ee , we do not have for the Coulombic correlation functional

E coul c [n] = E c [n; w coul
ee ] the same decomposition as in Eq. ( 4)

E coul c [n] = E lr,µ c [n] + E sr,µ c [n]. (5)
Instead, one can write E coul c

[n] exactly as

E coul c [n] = E lr,µ c [n] + E sr,µ c [n] + E lr-sr,µ c [n], (6)
defining the mixed long-range/short-range correlation functional E lr-sr,µ c

[n] which encompasses all the terms stemming from the non-linearity of E c [n; w ee ] with respect to w ee . By associating the mixed term E lr-sr,µ c

[n] to either the long-or the short-range part of the correlation energy, one obtains two possible long-range/shortrange decompositions. The first one is

E coul c [n] = E lr,µ c [n] + Ēsr,µ c [n], (7)
defining the new short-range correlation functional Ēsr,µ

c [n] = E sr,µ c [n] + E lr-sr,µ c [n]
. The second one is , we show now the expressions of these functionals in secondorder Görling-Levy perturbation theory [23]. We start from the second-order Coulombic correlation energy

E coul c [n] = Ēlr,µ c [n] + E sr,µ c [n], (8)
E coul,(2) c = i | Φ| Ŵ coul ee -V coul Hx |Φ i | 2 E s -E s,i , (9)
where V coul Hx = n(r)δE coul Hx /δn(r)dr is the Coulombic Hartree-exchange potential operator expressed with density operator n(r) and with the Coulombic Hartreeexchange functional E coul Hx [n], Φ and E s are the KS wave function and energy, Φ i and E s,i are the excited KS eigenfunctions and eigenvalues. Applying in Eq. ( 9) the longrange/short-range decomposition on Ŵ coul ee and V coul Hx and expanding leads to the second-order expressions for E lr,µ c , E sr,µ c and E lr-sr,µ c :

• The second-order long-range correlation energy writes

E lr,µ, (2) c
= i | Φ| Ŵ lr,µ ee -V lr,µ Hx |Φ i | 2 E s -E s,i , (10)
where V lr,µ Hx = n(r)δE lr,µ Hx /δn(r)dr is the longrange Hartree-exchange potential operator expressed with the long-range Hartree-exchange functional E lr,µ Hx [n] defined in the same way as the longrange exchange functional.

• Likewise, the second-order short-range correlation energy writes

E sr,µ, (2)
c = i | Φ| Ŵ sr,µ ee -V sr,µ Hx |Φ i | 2 E s -E s,i , (11)
where V sr,µ Hx = n(r)δE sr,µ Hx /δn(r)dr is the shortrange Hartree-exchange potential operator expressed with the short-range Hartree-exchange functional E sr,µ Hx [n] defined in the same way as the short-range exchange functional.

• Finally, the second-order mixed long-range/shortrange correlation functional is

E lr-sr,µ,(2) c = 2 i Φ| Ŵ lr,µ ee -V lr,µ Hx |Φ i Φ i | Ŵ sr,µ ee -V sr,µ Hx |Φ E s -E s,i . (12)

III. LOCAL DENSITY APPROXIMATION

For a given electron-electron interaction, the local density approximation to the previously introduced exchange and correlation functionals consists in locally transferring the corresponding energy of a uniform electron gas with the same interaction and with density equal to the local value of the inhomogeneous density.

Let us first consider the short-range functionals. The short-range exchange LDA functional associated to E sr,µ

x [n] writes

E sr,µ x,LDA [n] = n(r)ε sr,µ x,unif (n(r))dr, ( 13
)
where ε sr,µ x,unif (n) is the exchange energy per particle of a uniform electron gas with interaction w sr,µ ee [START_REF] Savin | Recent Developments of Modern Density Functional Theory[END_REF][START_REF] Toulouse | [END_REF]. Similarly, the short-range correlation LDA functional for Ēsr,µ

c [n] is Ēsr,µ c,LDA [n] = n(r)ε sr,µ c,unif (n(r))dr, ( 14
)
where εsr,µ c,unif (n) = ε coul c,unif (n) -ε lr,µ c,unif (n) is obtained as the difference of the Coulomb correlation energy, ε coul c,unif (n), and of the long-range correlation energy of a uniform electron gas with interaction w lr,µ ee [START_REF] Savin | Recent Developments of Modern Density Functional Theory[END_REF][START_REF] Toulouse | [END_REF]. The LDA for

E sr,µ c [n] writes E sr,µ c,LDA [n] = n(r)ε sr,µ c,unif (n(r))dr, (15)
where ε sr,µ c,unif is the correlation energy per particle of a uniform electron gas with interaction w sr,µ ee [26]. The LDA corresponding to the long-range functionals are obtained by difference to the Coulombic case:

E lr,µ x,LDA [n] = E coul x,LDA [n] -E sr,µ x,LDA [n], E lr,µ c,LDA [n] = E coul c,LDA [n] -Ēsr,µ c,LDA [n] and Ēlr,µ c,LDA [n] = E coul c,LDA [n] - E sr,µ c,LDA [n].

IV. RESULTS FOR THE HE ATOM

We consider the simple example of the He atom. We compare the accuracy of the LDA for the long-range or short-range energies as follows. For each value of mu (selected in [0, ∞[ ), we obtain an accurate value of the energy (see Refs. 17, 27 for details), we calculate the LDA error on the energy, and report in the plots the LDA error (ordinate) for each accurate energy (abscissa). We choose this way of plotting in order to emphasize the importance of the error for a given energy value to be recovered by the approximation: a good functional would yield small errors even for large contributions of the energy. Furthermore, this representation gives a unique scale for both long-and short-range expressions. The left end of the plots corresponds to vanishing functionals (i.e., µ = 0 for long-range functionals and µ → ∞ for short-range functionals). The right end of the plots corresponds to the usual (Coulombic) LDA case.

Let us consider first the LDA errors on E sr,µ

x and E lr,µ x , represented in Fig. 1. For all values of the exchange energy (except of course at the end points), the LDA error on E sr,µ x is systematically smaller than the error on E lr,µ x (the difference is of the order of 0.025 a.u. at an intermediate energy of -0.5 a.u.). Using the LDA for the short-range contribution to the exchange energy, rather than for the long-range contribution, allows, for a given energy, to do a smaller error on this energy, or equivalently, for a given error, to treat a larger part of the exchange energy. Note, in particular, that the LDA error on E sr,µ x becomes vanishingly small toward the left end of the plot, corresponding to a very-short-range interaction. This is in agreement with the fact that the LDA becomes exact in the limit of a very short-range interaction [17,28,29]. and Ēlr,µ c . It clearly appears that, in the whole energy range (except at the end points), the LDA errors on the short-range contributions are much smaller (in absolute value) than those on the long-range contributions. This confirms the appropriateness of the LDA for short-range correlations, as often pointed out in the literature (see, e.g., Refs. 17,20,21). Note that, contrary to the exchange case, the differences in the relative errors on the long-range and short-range contributions are important. For example, at an intermediate value of the energy of -0.02 a.u., the LDA errors on E lr,µ c and Ēlr,µ c are -0.0697 a.u. and -0.0673 a.u., respectively, while the LDA errors on E sr,µ c and Ēsr,µ c are as small as -0.0017 a.u. and -0.0001 a.u., respectively. One sees in addition that the LDA error on Ēsr,µ c is always significantly smaller (in absolute value) that the error on E sr,µ c . We conclude that the functional Ēsr,µ c is the best suited for a local density approximation.

V. CONCLUSIONS

As regards the long-range/short-range decomposition of the exchange functional of Eq. ( 4), the results of this work suggest that the LDA is more accurate for the shortrange contribution rather than for the long-range contribution.

Concerning the long-range/short-range decomposition of the correlation functional according to Eq. ( 7) or to Eq. ( 8), this work confirms that the LDA is more accurate for the short-range contributions. The presented results suggest in addition that the short-range correlation functional appearing in the decomposition of Eq. ( 7) is better suited for the LDA that the other short-range correlation functional appearing in the decomposition of Eq. ( 8).

In the context of the long-range/short-range decomposition in DFT, the present paper gives clues about the part of the energy, better adapted to (local) density functional approximations. We hope that it will incite to more systematic studies of this topic.

  defining the new long-range correlation functional Ēlr,µ c [n] = E lr,µ c [n] + E lr-sr,µ c [n]. All the long-range functionals vanish for µ = 0, E lr,µ=0 x = E lr,µ=0 c = Ēlr,µ=0 c = 0, and reduce to the Coulombic functionals for µ → ∞, E lr,µ→∞ x = E coul x and E lr,µ→∞ c = Ēlr,µ→∞ c = E coul c . Symmetrically, all the short-range functionals reduce to the Coulombic functionals for µ = 0, E sr,µ=0 x for µ → ∞, E sr,µ→∞ x = E sr,µ→∞ c = Ēsr,µ→∞ c = 0. The mixed long-range/shortrange correlation functional vanishes for both µ = 0 and µ → ∞, E lr-sr,µ=0 c = E lr-sr,µ→∞ c = 0. At first sight, the physical meaning of the decomposition of Eq. (6) may not be obvious. To get more insight on the nature of the terms E lr,µ c , E sr,µ c and E lr-sr,µ c

FIG. 1 :

 1 FIG. 1: LDA errors on E sr,µ x (solid curve) and E lr,µ x (dashed curve) with respect to -E sr,µ x and -E lr,µ x respectively, for the He atom.

Fig. 2

 2 compares the LDA errors on Ēsr,µ c , E sr,µ c , E lr,µ c

FIG. 2 :

 2 FIG. 2: LDA errors on Ēsr,µ c (thin solid curve), E sr,µ c (thick solid curve), E lr,µ c (thin dashed curve) and Ēlr,µ c (thick dashed curve) with respect to -Ēsr,µ c , -E sr,µ c , -E lr,µ c and -Ēlr,µ c

Acknowledgments

It is our pleasure to dedicate this paper to Annick Goursot, for courageously defending DFT over the years, building the environment which allowed us to pursue the present work.

We would also like to thank G. E. Scuseria for carefully reading our manuscript.


Julien Toulouse
email: toulouse@lct.jussieu.fr

Paola Gori-Giorgi

Andreas Savin
email: savin@lct.jussieu.fr

Scaling relations, virial theorem and energy densities for long-range and short-range density functionals

We analyze a decomposition of the Coulomb electron-electron interaction into a long-range and a short-range part in the framework of density functional theory, deriving some scaling relations and the corresponding virial theorem. We study the behavior of the local density approximation in the high-density limit for the long-range and the short-range functionals by carrying out a detailed analysis of the correlation energy of a uniform electron gas interacting via a long-range only electronelectron repulsion. Possible definitions of exchange and correlation energy densities are discussed and clarified with some examples.

I. INTRODUCTION

In recent years, there has been a growing interest in approaches of density functional theory (DFT) [1] based on a long-range/short-range decomposition of the Coulomb electron-electron interaction (see, e.g., Refs. [2][3][4][5][6]. The idea is to use different, appropriate approximations for the long-range and the short-range contributions to the usual density functionals of the Kohn-Sham (KS) scheme [7].

In this paper, we consider one possible longrange/short-range decomposition of the Coulombic density functionals. To gain insight into the introduced longrange and short-range density functionals, we examine some scaling relations, the virial theorem and energy densities. The local density approximation (LDA) appropriately defined for these functionals is also analyzed.

The long-range and short-range density functionals considered in this work are defined as follows (see also Refs. 5,[START_REF] Savin | Recent Developments of Modern Density Functional Theory[END_REF][START_REF] Leininger | [END_REF][10][11][12][13]. The starting point is the decomposition of the electron-electron Coulomb interaction w ee (r) = 1/r as w ee (r) = w lr,µ ee (r) + w sr,µ ee (r), (1) with the long-range interaction w lr,µ ee (r) = erf(µr)/r and the complement short-range part w sr,µ ee (r) = 1/rerf(µr)/r. The parameter µ controls the range of the decomposition. For µ = 0 the long-range interaction vanishes, w lr,µ=0 ee (r) = 0, while for µ → ∞ it reduces to the Coulomb interaction, w lr,µ→∞ ee (r) = 1/r. The Coulombic universal density functional F [n] is consequently decomposed as

F [n] = F lr,µ [n] + F sr,µ [n], (2)
where F lr,µ [n] is a long-range universal density functional and F sr,µ [n] is its short-range complement. The long-range functional is defined in the constrained-search formulation [14] by

F lr,µ [n] = min Ψ→n Ψ| T + Ŵ lr,µ ee |Ψ , ( 3
)
where T is the kinetic energy operator and Ŵ lr,µ ee = (1/2) n2 (r 1 , r 2 )w lr,µ ee (r 12 )dr 1 dr 2 the long-range interaction operator, expressed with the pair-density operator n2 (r 1 , r 2 ). The minimizing (multi-determinantal) wave function in Eq. ( 3) is denoted by Ψ lr,µ [n]. Introducing the non-interacting kinetic energy functional,

T s [n] = Φ[n]| T |Φ[n] , where Φ[n] is the KS determinant, F lr,µ [n] is written as F lr,µ [n] = T s [n] + E lr,µ Hxc [n], (4)
where the long-range Hartree-exchange-correlation functional, E lr,µ Hxc

[n] = E lr,µ H [n] + E lr,µ x [n] + E lr,µ c [n],
is the sum of the long-range Hartree functional

E lr,µ H [n] = 1 2 n(r 1 )n(r 2 )w lr,µ ee (r 12 )dr 1 dr 2 , (5)
the long-range exchange functional

E lr,µ x [n] = Φ[n]| Ŵ lr,µ ee |Φ[n] -E lr,µ H [n], (6)
and the long-range correlation functional

E lr,µ c [n] = Ψ lr,µ [n]| T + Ŵ lr,µ ee |Ψ lr,µ [n] -Φ[n]| T + Ŵ lr,µ ee |Φ[n] . (7)
The short-range functional, F sr,µ [n] = Ēsr,µ

Hxc [n] = E sr,µ H [n]+E sr,µ x [n]+ Ēsr,µ c [n],
is the sum of the short-range Hartree, exchange and correlation functionals, defined by complementarity to the Coulombic Hartree, exchange and correlation functionals,

E H [n], E x [n] and E c [n], E sr,µ H [n] = E H [n] -E lr,µ H [n], (8)
E sr,µ x [n] = E x [n] -E lr,µ x [n], (9)
Ēsr,µ

c [n] = E c [n] -E lr,µ c [n]. (10)
The LDA approximation to the long-range exchangecorrelation functional E lr,µ xc,LDA [n] is constructed from the exchange-correlation energy per electron ε lr,µ xc,unif of a uniform electron gas interacting with potential w lr,µ ee

E lr,µ xc,LDA [n] = n(r) ε lr,µ xc,unif (n(r)) dr. (11)
Similarly, the short-range LDA exchange-correlation functional Ēsr,µ xc,LDA [n] is defined with the complementary quantity ε sr,µ xc,unif = ε xc,unif -ε lr,µ xc,unif (see Refs. 8,15). The paper is organized as follows. In Sec. II, we examine some scaling properties of the long-range and shortrange functionals, and we discuss the high-density limit of the correlation functionals in LDA. In Sec. III, we derive the virial theorem satisfied by the long-range and short-range functionals. In Sec. IV, we examine longrange and short-range energy densities. Sec. V contains our conclusions.

All the relations derived in this work are more generally true for an interaction of the form w lr,µ ee (r) = β(µr)/r where β is a continuous and differentiable function satisfying β(x → ∞) = 1.

Atomic units (a.u.) are used throughout this work.

II. SCALING RELATIONS

In this section, we generalize some usual scaling relations of the Coulombic density functionals [16,17] (see also Ref. 18) to the case of the long-range and shortrange density functionals. The scaled wave function of a N -electron system corresponding to a uniform scaling of the electron coordinates by the scale factor γ > 0 is defined by (see, e.g., Refs. 16,[19][20][START_REF] Parr | Density-Functional Theory of Atoms and Molecules[END_REF] Ψ γ (r 1 , ..., r N ) = γ 3N/2 Ψ(γr 1 , ..., γr N ).

(

) 12
The wave function Ψ γ yields the scaled density n γ n γ (r) = γ 3 n(γr).

A. Scaling relation for the Hartree and exchange functionals

It can be easily verified from Eq. ( 5) that the longrange Hartree functional satisfies the following scaling relation

E lr,µγ H [n γ ] = γE lr,µ H [n]. ( 14
)
The KS determinant associated to the scaled density n γ is identical to the KS determinant associated to the density n after uniform scaling of the coordinates

Φ[n γ ] = Φ γ [n], (15)
and thus Φ

[n γ ]| Ŵ lr,µγ ee |Φ[n γ ] = Φ[n]| Ŵ lr,µ ee |Φ[n]
, leading to the same scaling relation for the long-range exchange functional

E lr,µγ x [n γ ] = γE lr,µ x [n]. (16)
The short-range Hartree and exchange functionals satisfy the same scaling relations

E sr,µγ H [n γ ] = γE sr,µ H [n], (17)
E sr,µγ x [n γ ] = γE sr,µ x [n]. (18)

B. Scaling relation for the correlation functionals

Let's introduce the following universal functional

F lr,µ,λ [n] = min Ψ→n Ψ| T + λ Ŵ lr,µ ee |Ψ , (19)
and the minimizing wave function is denoted by Ψ lr,µ,λ [n]. The scaled wave function Ψ lr,µ,λ γ

[n] gives the density n γ and minimizes Ψ|(

1/γ 2 ) T + (λ/γ) Ŵ lr,µγ ee |Ψ since Ψ lr,µ,λ γ |(1/γ 2 ) T + (λ/γ) Ŵ lr,µγ ee |Ψ lr,µ,λ γ = Ψ lr,µ,λ | T +λ Ŵ lr,µ
ee |Ψ lr,µ,λ . Consequently, doing the substitution λ → λ/γ, the wave function Ψ lr,µ,λ/γ γ

[n] gives the density n γ and minimizes (1/γ 2 ) Ψ| T + λ Ŵ lr,µγ ee |Ψ . Therefore, we conclude that

Ψ lr,µγ,λ [n γ ] = Ψ lr,µ,λ/γ γ [n]. (20)
Applying this relation with λ = 1, we find the following scaling relation for the long-range universal functional

F lr,µγ [n γ ] = γ 2 F lr,µ,1/γ [n], (21)
and consequently for the long-range correlation functional

E lr,µγ c [n γ ] = γ 2 E lr,µ,1/γ c [n]. ( 22
)
where

E lr,µ,λ c [n] = Ψ lr,µ,λ [n]| T + λ Ŵ lr,µ ee |Ψ lr,µ,λ [n] - Φ[n]| T +λ Ŵ lr,µ ee |Φ[n]
. The short-range correlation functional satisfies the same scaling relation

Ēsr,µγ c [n γ ] = γ 2 Ēsr,µ,1/γ c [n], (23)
where Ēsr,µ,λ

c [n] = E λ c [n] -E lr,µ,λ c [n] and E λ c [n]
is the Coulombic correlation functional along the linear adiabatic connection. To our knowledge, Eq. ( 20) is new while Eq. ( 21) has already been mentioned by Yang [START_REF] Yang | [END_REF].

C. High-density limit of the correlation functionals

We study now the long-range and short-range correlation functionals under uniform scaling to the high-density limit (γ → ∞). If the ground-state of the KS system

is non-degenerate, E lr,µγ c [n γ ] goes to a constant when γ → ∞ lim γ→∞ E lr,µγ c [n γ ] = Ēlr,µ,(2) c [n], (24)
where E lr,µ,(2) c

[n] = (1/2)(∂ 2 E lr,µ,λ c

[n]/∂λ 2 ) λ=0 is the second-order correlation energy in the Görling-Levy perturbation theory [23], just as in the Coulombic case [17,24]. We have a similar behavior for Ēsr,µγ

c [n γ ] lim γ→∞ Ēsr,µγ c [n γ ] = Ēsr,µ,(2) c [n], (25)
with Ēsr,µ,( 2)

c [n] = (1/2)(∂ 2 Ēsr,µ,λ c [n]/∂λ 2 ) λ=0 .
It is interesting to study how the long-range and shortrange LDA functionals behave in the high-density limit, and to understand whether they can overcome some of the well-known problems of standard LDA in this regime (see, e.g., [25]). For the uniform electron gas of density n, the scaling transformation of Eq. ( 13) is simply written as r s → r s /γ, where r s = (4πn/3) -1/3 is the Wigner-Seitz radius. The correlation energy per particle of the Coulombic electron gas diverges in the high-density limit logarithmically [26]

ε c,unif (r s /γ) ∼ γ→∞ -A ln γ + • • • , (26)
where A = (1 -ln 2)/π 2 . The Coulombic LDA functional thus cannot recover the correct high-density scaling for systems with a non-degenerate ground state [25]. The origin of the divergence of Eq. ( 26) has been widely analyzed in terms of the electron-gas pair density, both in real and reciprocal space (see, e.g., [27][28][29]), and, more recently, in terms of the momentum distribution [30].

The investigation of the γ → ∞ performances of the long-and short-range LDA functionals implies the study of the high-density behavior of the correlation energy of an electron gas interacting with w lr,µ ee . Following Refs. 27, 28, we consider for this "long-range" electron gas the Fourier transform of the pair density, the static structure factor S lr,µ (k, r s ). The correlation energy per particle, ε lr,µ c,unif , is given in terms of this function as

ε lr,µ c,unif (r s ) = 1 2 1 0 dλ dk (2π) 3 wlr,µ ee (k)S lr,µ,λ c (k, r s ), (27)
where S lr,µ,λ c (k, r s ) = S lr,µ,λ (k, r s )-S KS (k, r s ) is the correlation part of the static structure factor for the system of density parameter r s interacting with λw lr,µ ee , and

wlr,µ ee (k) = 4 π k 2 e -k 2 /(4 µ 2 ) ≡ 4 π k 2 η(k/µ) (28)
is the Fourier transform of the long-range interaction w lr,µ ee (r) = erf(µr)/r. The results derived below hold more generally for any long-range interaction whose Fourier transform is of the form (28), 4 π k 2 η(k/µ), with η(x → 0) = 1, i.e., for any long-range interaction whose real space form is β(µr)/r with β(x → ∞) = 1, as stated at the end of the Introduction.

For the Coulomb gas, the random-phase approximation (RPA) provides an expression for S λ c which is exact for small (≪ 1) values of the scaled variable q = k/2k F = αr s k/2 (with α 3 = 4/9π) and which gives the exact correlation energy when r s → 0 (see, e.g., Refs. 28,31). The arguments for the validity of RPA in the high-density limit [31] can be extended to an interaction of the kind (28) as long as, when r s → 0, wlr,µ ee (k) diverges for small k as k -2 : in this case, in fact, the perturbation series expansion for the correlation energy contains as leading term an infinite number of divergent direct diagrams whose resummation yields the RPA expression for the correlation energy, as in the Coulomb gas. The RPA S lr,µ,λ c reads S lr,µ,λ c,RPA (q, r s ) = -

6q k F η 2q α rs µ × ∞ 0 λ χ0 (q, iu) 2 du q 2 -λ α r s χ0 (q, iu) η 2q α rs µ , (29)
where χ0 (q, iu) = (π/k F )χ 0 (q, iu, r s ) is a reduced noninteracting response function, expressed in terms of the reduced imaginary frequency u = -iω/(2qk 2 F ), that does not contain any explicit dependence on r s . In the case of the Coulombic gas (which corresponds to µ = ∞, that is, η = 1), if one simply lets r s go to zero, the r.h.s. of Eq. ( 29) gives a static structure factor which behaves like 1/q for small q, thus yielding the logarithmic divergence of the correlation energy. This divergence comes from the combination of the two limits, r s → 0 and q → 0, and requires an appropriate treatment, for which it is convenient to divide the integral of Eq. ( 27) in two parts, q1 0 dq + ∞ q1 dq, with q 1 ≪ 1. The second part of this integral is finite when r s → 0, while the first part yields the logarithmic divergence and it is better studied by switching to the scaled variable y = q/ √ r s [27,28,30], which gives an exact scaling for S lr,µ,λ c when q ≪ 1 and r s → 0,

S lr,µ,λ c (q ≪ 1, r s → 0) = -6 √ r s α y η( 2 √ rs y α µ rs ) ∞ 0 λ R(u) 2 du y 2 -λ α R(u) η( 2 √ rs y α µ rs ) , (30)
where R(u) = (u arctan 1 u -1)/π. Equation (30) shows that if we rescale µ, µ → µ/r s , the factors η( 2 √ rs y α µ rs ) become equal to 1 when r s = 0, so that the small-q part of S lr,µ/rs,λ c scales exactly to the same limit of the Coulombic gas, i.e., S lr,µ/rs,λ

c (q ≪ 1, r s ) ∼ rs→0 -6 √ r s α f (y, λ), (31)
where the function f (y, λ) [27,28] behaves as y for small y and as 1/y for large y, and is reported in Fig. 1 in the case λ = 1. Then everything goes as in the Coulomb gas [28]: in the small-q part of the integral in Eq. ( 27) we can replace S lr,µ/rs,λ c

with Eq. ( 31), obtaining an expression of the kind q1/ √ rs 0 dyf (y, λ). When r s → 0, even if q 1 is small, the upper limit of this integral is large, and the 1/y behavior of f (y ≫ 1) causes the logarithmic divergence of the correlation energy per particle as γ → ∞

ε lr,µγ c,unif (r s /γ) ∼ γ→∞ -A ln γ + • • • , (32)
with exactly the same A of Eq. ( 26). The short-range correlation energy per particle of the uniform electron gas is just the difference between the correlation energies of the Coulomb gas and of the "longrange" gas, so that an exact cancellation of the logarithmic term occurs if we rescale µ,

lim γ→∞ εsr,µγ c,unif (r s /γ) = h(µ, r s ). ( 33
)
The function h(µ, r s ) remains finite as long as µr s is greater than zero. This means that one can improve the LDA performances in the high-density limit by rescaling µ either locally, µ ∝ 1/r s (r), or globally by choosing a µ bigger than 1/r min s , where r min s is the minimum value of r s (r) in the given system. In other words, the rescaled short-range functional allows to perform a cut-off of the correlation transferred from an electron gas, i.e., we can eliminate the long-range correlations that occur in an electron gas but that do not occur in a confined system.

III. VIRIAL THEOREM

In this section, we generalize the virial theorem of the Coulombic density functionals [16] (see also Refs. 20,32) to the case of the long-range and short-range density functionals.

The wave functions of the type Ψ lr,µ γ [n 1/γ ] give the density n, independently of γ, and can therefore be used as trial wave functions in the variational definition of F lr,µ [n] [Eq. ( 3)]. As Ψ lr,µ γ [n 1/γ ] reduces to the minimizing wave function Ψ lr,µ [n] at γ = 1, the stationarity condition implies that the derivative with respect to γ vanishes at γ = 1

d dγ Ψ lr,µ γ [n 1/γ ]| T + Ŵ lr,µ ee |Ψ lr,µ γ [n 1/γ ] γ=1 = 0. ( 34
)
The kinetic and electron-electron interaction energies have simple uniform coordinate scalings

d dγ (γ 2 T lr,µ [n 1/γ ] + γW lr,µ/γ ee [n 1/γ ]) γ=1 = 0, ( 35
)
where

T lr,µ [n] = Ψ lr,µ [n]| T |Ψ lr,µ [n] and W lr,µ ee [n] = Ψ lr,µ [n]| Ŵ lr,µ ee |Ψ lr,µ [n] .
Performing the derivative with respect to γ leads to

2T lr,µ [n] + W lr,µ ee [n] -µ ∂W lr,µ ee [n] ∂µ = d dγ (T lr,µ [n γ ] + W lr,µ ee [n γ ]) γ=1 . ( 36
)
Using the virial relation for non-interacting kinetic functional [18] 2T

s [n] = d dγ T s [n γ ] γ=1 , (37)
Eq. ( 36) simplifies to

T lr,µ c [n] + E lr,µ Hxc [n] -µ ∂E lr,µ Hxc [n] ∂µ = d dγ E lr,µ Hxc [n γ ] γ=1 = -n(r)r.∇ δE lr,µ Hxc [n] δn(r) dr, (38)
where T lr,µ c

[n] = T lr,µ [n] -T s [n] and the last equality has been obtained through a integration by parts.

It can be verified that the long-range Hartree, exchange and correlation functionals obey separate virial relations, just as for the Coulombic case,

E lr,µ H [n] -µ ∂E lr,µ H [n] ∂µ = -n(r)r.∇ δE lr,µ H [n] δn(r) dr, (39)
E lr,µ x [n] -µ ∂E lr,µ x [n] ∂µ = -n(r)r.∇ δE lr,µ x [n] δn(r) dr, (40)
T lr,µ c [n]+E lr,µ c [n]-µ ∂E lr,µ c [n] ∂µ = -n(r)r.∇ δE lr,µ c [n] δn(r)
dr.

(41) The same relations are also valid for the short-range Hartree, exchange and correlation functionals

E sr,µ H [n] -µ ∂E sr,µ H [n] ∂µ = -n(r)r.∇ δE sr,µ H [n] δn(r) dr, (42)
E sr,µ x [n] -µ ∂E sr,µ x [n] ∂µ = -n(r)r.∇ δE sr,µ x [n] δn(r) dr, (43)
T sr,µ

c [n]+ Ēsr,µ c [n]-µ ∂ Ēsr,µ c [n] ∂µ = -n(r)r.∇ δ Ēsr,µ c [n] δn(r) dr, (44) where
T sr,µ c [n] = T c [n] -T lr,µ c [n] and T c [n] is the usual Coulombic correlation kinetic functional.
For the corresponding virial relations in the uniform electron gas, see Ref. 33.

IV. ENERGY DENSITIES

In this section, we examine long-range and short-range energy densities defined from pair densities or from the virial theorem. Energy densities are always useful to analyze approximations or derive new approximations (see, e.g., Refs. [34][35][36][37][38][39][40].

A. Energy densities from pair densities

In Ref. 41, energy densities for modified interactions defined from pair densities have been discussed. We recall that an energy density associated to the long-range exchange energy E lr,µ x can be defined by e lr,µ

x,pd (r) =

1 2 n 2,x (r, r 12 )w lr,µ ee (r 12 )dr 12 , (45)
where n 2,x (r, r 12 ) is the KS exchange pair density. Likewise, the corresponding short-range energy density writes

e sr,µ x,pd (r) = 1 2 n 2,x (r, r 12 )w sr,µ ee (r 12 )dr 12 . (46)
The definition of correlation energy densities from pair densities requires an integration over an adiabatic connection. For instance, an energy density associated to the long-range correlation functional E lr,µ c can be written as e lr,µ c,pd (r) = These energy densities involves two-electron quantities that can complicate their evaluation.

B. Energy densities from the virial theorem

Long-range and short-range energy densities can be defined from the virial theorem, just as for the Coulombic case (see, e.g., Refs. 34,38). The virial relation of Eq. ( 40) leads indeed to the following long-range exchange energy density e lr,µ

x,vir (r) = -µ

∞ µ dξ ξ 2 n(r) r • ∇v lr,ξ x (r), (49)
where v lr,µ x (r) = δE lr,µ x [n]/δn(r). Likewise, Eq. ( 43) leads to the following short-range exchange energy density e sr,µ

x,vir (r) = -µ

∞ µ dξ ξ 2 n(r) r • ∇v sr,ξ x (r), (50)
where v sr,µ x (r) = δE sr,µ x [n]/δn(r). The virial relation of Eq. ( 41

where v sr,µλ,λ c (r) = δ Ēsr,µλ,λ c

[n]/δn(r). These energy densities have the advantage of involving only one-electron quantities.

C. Results on the He atom

As a simple illustration, we have calculated for the He atom the long-range and short-range energy densities e lr,µ

x,pd (r), e sr,µ x,pd (r), e lr,µ x,vir (r) and e sr,µ x,vir (r) for µ = 1 a.u., as well as the Coulombic energies densities e x,pd (r) = e lr,µ

x,pd (r) + e sr,µ x,pd (r) and e x,vir (r) = e lr,µ x,vir (r) + e sr,µ x,vir (r). For a two-electron system, the exchange pair density and exchange potentials are directly deducible from the density. Using an accurate density, accurate exchange energy densities are thus easily obtained.

Fig. 2 compares the accurate and LDA radial exchange energy densities as a function of the distance to the nucleus r. One sees that the energy densities defined from the exchange pair density and from the virial theorem are qualitatively similar. At small r (r 0.5 a.u.), the LDA slightly overestimates the accurate Coulombic energy densities. At large r (r 0.5 a.u.), the LDA importantly underestimates the accurate Coulombic energy densities. The contribution at large r remains important in the long-range energy densities while it is significantly reduced in the short-range energy densities. Consequently, the LDA better performs for the short-range energy densities. x,pd (r) [Eq. ( 45)] and 4πr 2 e lr,µ x,vir (r) [Eq. (49)], (c) short-range radial energy densities 4πr 2 e sr,µ

x,pd (r) [Eq. ( 46)] and 4πr 2 e sr,µ

x,vir (r) [Eq. ( 50)], with µ = 1 a.u..

In the case of the energy density defined from the exchange pair density, the better performance of the LDA at small distance r can be easily explained in term of the spherical average of the exchange hole n x (r, r 12 ) = n 2,x (r, r 12 )/n(r) n sph. avr.

x (r, r 12 ) = 1 4π n x (r, r 12 )dΩ r12 ,

which is represented for the He atom in Fig. 3 with re-spect to the interelectronic distance r 12 for two positions of the reference electron r = 0 and r = 1 a.u.. For r = 0, both the accurate and LDA exchange holes are centered at r 12 = 0, making the LDA a reasonable approximation. For r = 1 a.u., the accurate hole is centered near r 12 = 1 a.u. while the LDA hole is still centered at r 12 = 0, leading to an important underestimation of the hole. x (r, r12) [Eq. ( 53)] for the He atom with r = 0 and r = 1 a.u..

V. CONCLUSIONS

We have analyzed a short-range and long-range decomposition of the Coulomb electron-electron interaction and we have derived some exact scaling relations for the corresponding density functionals. The study of the LDA approximation has shown that in the high-density limit the short-range functional scales to a constant, thus opening the possibility of ameliorating the performances of the Coulomb LDA functional in this regime. Possible definitions of energy densities obtained from pair densities and from the virial theorem have been presented. Results on the He atom suggest that the LDA approximation can give accurate short-range exchange energy densities.

FIG. 1 :

 1 FIG.1:The function f (y, λ) of Eq. (31) at λ = 1.

  (r, r 12 ) ∂w lr,ξ ee (r 12 ) ∂ξ dr 12 , (47) where n lr,ξ 2,c (r, r 12 ) is the correlation pair density for the long-range interaction w lr,ξ ee (r 12 ). The corresponding short-range correlation energy density is ēsr,µ c,pd (r(r, r 12 ) ∂w lr,ξ ee (r 12 ) ∂ξ dr 12 . (48)

1 dλ λ 3 1 dλ λ 3

 1313 ) generalized to the linear adiabatic connection, T lr,µc [n] + E lr,µ c [n]µ ∂E lr,µ c [n]/∂µ = -n(r)r.∇δE lr,µ c [n]/δn(r)dr where T lr,µ c [n] is given by T lr,µ,λ c [n] = E lr,µ,λ c [n]λ ∂E lr,µ,λ c [n]/∂λ (see Refs.38, 42) enables to define a long-range correlation energy density e lr,µ c,vir (r) = -∞ n(r) r • ∇v lr,µλ,λ c (r), (51) where v lr,µλ,λ c (r) = δE lr,µλ,λ c [n]/δn(r). Likewise, the virial relation of Eq. (44) leads to the short-range correlation energy density ēsr,µ c,vir (r) = -∞ n(r) r • ∇v sr,µλ,λ c (r),

FIG. 2 :

 2 FIG.2: Accurate (solid curves) and LDA (dashed curves) radial exchange energy densities defined from pair densities (thick curves) and from the virial theorem (thin curves) for the He atom: (a) Coulombic radial energy densities 4πr 2 e x,pd (r) and 4πr 2 ex,vir(r) (b) long-range radial energy densities 4πr 2 e lr,µx,pd (r) [Eq. (45)] and 4πr 2 e lr,µ x,vir (r) [Eq. (49)], (c) short-range radial energy densities 4πr 2 e sr,µx,pd (r) [Eq. (46)] and 4πr 2 e sr,µx,vir (r) [Eq. (50)], with µ = 1 a.u..

FIG. 3 :

 3 FIG.3: Accurate (solid curves) and LDA (dashed curves) spherical-average exchange hole n sph. avr.


Julien Toulouse

Paola Gori-Giorgi

Andreas Savin

A short-range correlation energy density functional with multi-determinantal reference

We introduce a short-range correlation density functional defined with respect to a multideterminantal reference which is meant to be used in a multi-determinantal extension of the Kohn-Sham scheme of density functional theory based on a long-range/short-range decomposition of the Coulomb electron-electron interaction. We construct the local density approximation for this functional and discuss its performance on the He atom. c and ∆ Ēsr,µ c,md = Ēsr,µ c,md,LDA -Ēsr,µ c,md with respect to µ for the He atom.

I. INTRODUCTION

One of the main difficulties in the Kohn-Sham (KS) [1] scheme of density functional theory (DFT) [2] is to find approximations for the exchange-correlation energy functional that correctly describe (near-)degeneracy or longrange (e.g., van der Waals) correlation effects. To circumvent this difficulty, a multi-determinantal extension of the KS scheme based on a long-range/short-range decomposition of the Coulomb electron-electron interaction has been proposed [3][START_REF] Savin | Recent Advances in Density Functional Theory[END_REF][START_REF] Savin | Recent Developments of Modern Density Functional Theory[END_REF][START_REF] Leininger | [END_REF][7][8][9][10][START_REF] Stoll | Density Functional Method in Physics[END_REF][START_REF] Toulouse | [END_REF]. The idea behind this separation is that correlation effects due to the short-range part, involving the correlation cusp, could well be described by the local density approximation (appropriately modified); correlation connected with the long-range part could well be dealt with using standard wave-function methods of quantum chemistry.

In this approach, the ground-state energy of a Nelectron system in a nuclei-electron potential v ne (r) is obtained in principle exactly by minimization over multideterminantal wave functions Ψ E = min

Ψ Ψ| T + Vne + Ŵ lr,µ ee |Ψ +E sr,µ H [n Ψ ] + E sr,µ x [n Ψ ] + Ēsr,µ c [n Ψ ] , (1)
where T is the kinetic energy operator, Vne = i v ne (r i ) is the nuclei-electron interaction operator, Ŵ lr,µ ee = i<j w lr,µ ee (r ij ) is a long-range electronelectron interaction operator with w lr,µ ee (r) = erf(µr)/r,

E sr,µ H [n] = (1/2) n(r 1 )n(r 2 )w sr,µ ee (r 12 )dr 1 dr 2 is a short- range Hartree functional with w sr,µ ee (r) = erfc(µr)/r, E sr,µ x [n] = Φ[n]| Ŵ sr,µ ee |Φ[n] -E sr,µ H [n] is a short-range exchange functional where Ŵ sr,µ ee = i<j w sr,µ ee (r ij ) and Φ[n] is the KS determinant, Ēsr,µ c
[n] is a short-range correlation functional defined so that Eq. ( 1) is exact, and n Ψ is the density coming from Ψ. The minimizing wave function in Eq. (1) will be denoted by Ψ lr,µ . In these * Electronic address: savin@lct.jussieu.fr equations, µ is a parameter controlling the range of the decomposition of the Coulomb interaction. In practice, approximations must be used for the wave function Ψ lr,µ and the short-range functionals E sr,µ

x [n] and Ēsr,µ c

[n]. In particular, the local density approximations have been constructed for E sr,µ

x [n] and Ēsr,µ c

[n] [START_REF] Savin | Recent Developments of Modern Density Functional Theory[END_REF]13]. For µ = 0, Eq. ( 1) reduces to the KS scheme. In fact, in this case, the long-range interaction vanishes, Ŵ lr,µ=0 ee = 0, and the short-range functionals E sr,µ=0 H

[n], E sr,µ=0

x

[n] and Ēsr,µ=0 [n] reduce to the Hartree, exchange and correlation functionals of the KS theory.

Former experience with Eq. ( 1) has shown that in general the quality of the wave-function Ψ lr,µ obtained with a given approximate functional is much better than that of the functional itself. To extract the maximum information from Ψ lr,µ , we propose in this work to compute the ground-state energy as

E = Ψ lr,µ | T + Vne + Ŵee |Ψ lr,µ + Ēsr,µ c,md [n Ψ lr,µ ], (2)
where Ŵee = i<j 1/r ij is the full Coulomb interaction operator and Ēsr,µ c,md [n] is a new short-range correlation functional defined such as Eq. ( 2) used with the exact wave function Ψ lr,µ is exact.

We note that Eq. (2) can formally be made selfconsistent by generalizing the "optimized effective potential" (OEP) approach (see, e.g., Refs. [START_REF]A Primer in Density Functional Theory[END_REF]15) to the multi-determinantal extension of the KS scheme

E = inf v Ψ lr,µ [v]| T + Vne + Ŵee |Ψ lr,µ [v] + Ēsr,µ c,md [n v ] , (3)
where the infinimum is search over one-electron potentials v(r), and Ψ lr,µ [v] and n v are, respectively, the ground-state multi-determinantal wave function and density of T + Ŵ lr,µ ee

+ i v(r i ). If it exists, the mini- mizing potential is v(r) = v ne (r) + δE sr,µ H [n]/δn(r) + δE sr,µ x [n]/δn(r) + δ Ēsr,µ c
[n]/δn(r), establishing the link with Eq. (1). Eq. ( 3) reduces to the OEP formulation of the KS theory when µ = 0. In this case, in fact, Ψ lr,µ=0 is the KS determinant and the functional Ēsr,µ=0 c,md [n] reduces to the correlation functional of the KS theory. In practice, because the potential corresponding to the functional derivative of Ēsr,µ c,md [n] is small, we expect that the corrections on Ψ lr,µ brought by self-consistency are negligible with respect to the errors due to the approximations on Ēsr,µ c,md [n] [16]. This work is devoted to the study of the functional Ēsr,µ c,md [n]: we turn our attention entirely to the correlation energy, without combining it with an approximate functional for exchange. It is thus a pleasure to dedicate this paper to Professor Hermann Stoll who has been a pioneer in the study and application of correlation energy density functionals [17,18].

The paper is organized as follows. In Sec II, we discuss the short-range correlation functional Ēsr,µ c,md [n] and its relation to the functional Ēsr,µ c

[n]. In Sec. III, we construct a local density approximation for Ēsr,µ c,md [n]. In Sec. IV, we assess the accuracy of this approximation for the He atom. Sec. V contains our conclusions. Atomic units (a.u.) are used throughout this work.

II. THE SHORT-RANGE CORRELATION FUNCTIONAL Ēsr,µ c,md [n]
The short-range correlation functional Ēsr,µ c,md [n] in Eq. ( 2) is defined with respect to the multi-determinantal wave function Ψ lr,µ , in contrast to the short-range correlation functional Ēsr,µ c

[n] in Eq. ( 1) defined with respect to the one-determinant wave function Φ. It is easy to see that these two functionals are related to each other through

Ēsr,µ c,md [n] = Ēsr,µ c [n] + ∆ lr-sr,µ [n], (4)
where

∆ lr-sr,µ [n] = -Ψ lr,µ [n]| Ŵ sr,µ ee |Ψ lr,µ [n] -Φ[n]| Ŵ sr,µ ee |Φ[n] . (5)
The quantity ∆ lr-sr,µ vanishes for µ = 0 and µ → ∞.

It is interesting to study the behavior of Ēsr,µ c,md [n] in the limit of a very short-range interaction, i.e. when µ → ∞. In this limit, the short-range interaction behaves as [START_REF] Toulouse | [END_REF] w sr,µ ee (r) =

π µ 2 δ(r) + O 1 µ 3 , (6)
leading to the following asymptotic expansion of Ēsr,µ

c [8, 12, 19] Ēsr,µ→∞ c = π 2µ 2 n 2,c (r, r)dr + O 1 µ 3 , (7)
where n 2,c (r, r) is the correlation on-top pair density. The asymptotic expansion of ∆ lr-sr,µ as µ → ∞ is obtained similarly from its definition, Eq. ( 5), leading to

∆ lr-sr,µ→∞ = - π 2µ 2 n 2,c (r, r)dr + O 1 µ 3 . ( 8
)
The first terms in Eqs. ( 7) and ( 8) cancel, and therefore Ēsr,µ c,md decays at least as 1/µ 3 when µ → ∞.

III. LOCAL DENSITY APPROXIMATION

A local density approximation (LDA) can be constructed for Ēsr,µ c,md [n] Ēsr,µ c,md,LDA [n] = n(r)ε sr,µ c,md,unif (n(r))dr, (9) where the corresponding correlation energy per particle in the uniform electron gas εsr,µ c,md,unif (n) is given by

εsr,µ c,md,unif (n) = εsr,µ c,unif (n) + ∆ lr-sr,µ unif (n). (10)
In Eq. ( 10), εsr,µ c,unif (n) is the correlation energy per particle defining the LDA approximation for Ēsr,µ c

[n] (see Refs. [START_REF] Savin | Recent Developments of Modern Density Functional Theory[END_REF]13), and ∆ lr-sr,µ unif (n) is given by

∆ lr-sr,µ unif (n) = - n 2 ∞ 0 g lr,µ c,unif (r, n)w sr,µ ee (r)4πr 2 dr, (11)
where g lr,µ c,unif (r, n) is the correlation pair-distribution function of a uniform electron gas with long-range interaction w lr,µ ee (r) and density n. The correlation hole of this "long-range" electron gas is then given by n g lr,µ c,unif (r, n). Since an estimate of the energy εsr,µ c,unif (n) from coupledcluster calculations is available [START_REF] Savin | Recent Developments of Modern Density Functional Theory[END_REF]13], we only need to compute the term ∆ lr-sr,µ unif (n) to build the LDA functional of Eq. ( 9). In order to estimate ∆ lr-sr,µ unif (n) we proceed as follows. We first notice that w sr,µ ee (r) = erfc(µr)/r in the integrand of Eq. ( 11) only samples the part of g lr,µ c,unif (r, n) corresponding to r 1/µ. For the standard uniform electron gas (with full interaction 1/r) the "extended Overhauser model" [20] proved to be able to yield accurate results for g c,unif (r, n) in the short-range region defined by r ≤ r s , where r s = (4π n/3) -1/3 . We can thus use this simple model to calculate g lr,µ c,unif (r, n) and to produce an estimate for ∆ lr-sr,µ unif (n) that should be reliable for µ-values for which µr s 1.

The scattering equations of the "extended Overhauser model" are widely explained in Refs. 20, 21. Here we solved the same equations with the electron-electron interaction erf(µr)/r screened by a sphere of radius r s of uniform positive charge density n and attracting the electrons with the same modified interaction,

V eff (r, r s , µ) = erf(µr) r - |r ′ |≤rs n erf(µ|r ′ -r|) |r ′ -r| dr ′ . ( 12
)
This potential is reported in the Appendix of Ref. 22, where it has been used for two-electron atoms with very accurate results for the corresponding short-range correlation energy. V eff (r, r s , µ) is a screened potential that tends to the "Overhauser potential" [20,23] when µ → ∞, and which goes to zero when µ → 0. As in the original work of Overhauser [23], the idea behind Eq. ( 12) is that the radius of the screening "hole" is exactly equal to r s . c,unif (r, n) of a uniform electron gas of density 3/(4πr 3 s ) (here rs = 2) with long-range-only interaction erf(µr)/r obtained from the "extended Overhauser model" [20] [see Eq. ( 12)]. The correlation hole is given by n g lr,µ c,unif (r, n).

A sample of the pair-correlation functions g lr,µ c,unif (r, n) that we have obtained is reported in Fig. 1. As long as µ is not large we clearly see the absence of the cusp [g ′ c,unif (r = 0) is not zero for a system with interaction 1/r at small r, but it is zero for the erf(µr)/r interaction]. As expected, as µ increases the hole deepens, and for very large µ we see that the cusp starts to appear.

Some of the values of ∆ lr-sr,µ unif (n) for µr s 1 are reported in Table I. An estimate of ∆ lr-sr,µ unif (n) in the region not accessible with the extended Overhauser model, µr s 1, has been obtained by a simple interpolation between our data and zero, since, as explained in Sec. II, ∆ lr-sr,µ vanishes when µ → 0. In the opposite limit, µ → ∞, ∆ lr-sr,µ unif (n) behaves as in Eq. ( 8), which for a system of uniform density reads ∆ lr-sr,µ→∞ unif

(n) = - 3 g c,unif (0, n) 8 r 3 s µ 2 + O 1 µ 3 , (13)
where g c,unif (0, n) is the on-top value (r = 0) of the paircorrelation function of the Coulombic uniform electron gas of density n [20]. We found that the ∆ lr-sr,µ unif (n) computed with the extended Overhauser model accurately recover this limiting behavior.

For future applications, a more accurate LDA functional for ∆ lr-sr,µ (especially for µr s 1) will be available from quantum Monte Carlo calculations [24].

IV. RESULTS FOR THE HE ATOM

For the He atom, the short-range correlation energies Ēsr,µ c and Ēsr,µ c,md have been calculated with a precision of the order of 1 mH as follows. An accurate density is calculated at the full configuration interaction level with a large Gaussian basis set and the optimization of the potential in the Legendre transform formulation [25,26] of density functionals enables to compute accurately the |, meaning that in Eq. ( 5) Ψ lr,µ | Ŵ sr,µ ee |Ψ lr,µ < Φ| Ŵ sr,µ ee |Φ which seems natural for a repulsive interaction. Tab. II compare the LDA errors on Ēsr,µ c and Ēsr,µ c,md . One sees that the LDA errors for this two short-range correlation energies are of the same order of magnitude for all values of µ.

V. CONCLUSIONS

In this work, we have reexamined the multideterminantal extension of the KS scheme based on a long-range/short-range decomposition of the Coulomb electron-electron interaction. Contrary to previous works where the short-range correlation functional was defined with respect to the KS determinant, we have introduced a new short-range correlation functional defined with respect to the multi-determinantal wave function. We have constructed the local density approximation for this new functional. The example of the He atom suggests that the local density approximation is essentially as accurate as for the short-range correlation functional defined with respect to the KS determinant. We believe that this work paves the way to a multi-determinantal extension of the KS scheme using a correlation-only density functional. 10)- [START_REF] Stoll | Density Functional Method in Physics[END_REF]] computed from the extended Overhauser model [20]. µ 0.00 0.25 0.50 0.75 1.00 1.50 2.00 3.00 5.00 ∆ Ēsr,µ c -0.071 -0.053 -0.035 -0.022 -0.014 -0.005 -0.001 0.000 0.000 ∆ Ēsr,µ c,md -0.071 -0.048 -0.030 -0.019 -0.013 -0.006 -0.002 -0.002 0.000

  c

2 FIG. 1 :

 21 FIG.1:A sample of the pair-correlation functions g lr,µ c,unif (r, n) of a uniform electron gas of density 3/(4πr 3 s ) (here rs = 2) with long-range-only interaction erf(µr)/r obtained from the "extended Overhauser model"[20] [see Eq. (12)]. The correlation hole is given by n g lr,µ c,unif (r, n).

FIG. 2 :

 2 FIG.2: Short-range correlation energies Ēsr,µ c (thin curves) and Ēsr,µ c,md (thick curves) with respect to µ for the He atom. Accurate calculations (solid curves) are compared to the LDA approximation (dashed curves).

  000219 10 3 5.63 • 10 -5 15 3 2.52 • 10 -5 2 5 0.000335 3 5 0.000151 5 4 0.000102 10 4 2.59 • 10 -5 15 4 1.16 • 10 -5 2 6 0.000202 3 6 9.08 • 10 -5 5 5 5.52 • 10 -5 10 5 1.39 • 10 -5 15 5 6.19 • 10 -6

TABLE I :

 I A sample of the values of ∆ lr-sr,µ unif [see Eqs. (

TABLE II :

 II LDA errors ∆ Ēsr,µ

c = Ēsr,µ c,LDA -Ēsr,µ