A Polyakov
email: andrey.polyakov@inria.fr

J.-P Richard

Suboptimal LMI-based Solution of Minimum Time Control Problem

The problem of time optimal control design is considered for a chain of integrators. The suboptimal solution based on Implicit Lyapunov Function (ILF) method is presented in the form of continuous stabilizing feedback. The Semi-Definite Programming (SDP) problem with the constraints in the form of Linear Matrix Inequalities (LMI) is obtained for tuning the optimal parameters. The suboptimal solution is compared with the minimum discontinuous feedback on numerical example.

INTRODUCTION

The control applications frequently asks for algorithms, which provide terminations of transition processes in a minimal possible time. Such problem statements usually appear in robotic systems, aerospace applications, underwater/surface vehicles control systems, etc. Despite of long history of the optimal control theory [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF], [START_REF] Bellman | Dynamic Programming[END_REF]) time optimal control problems are still subjects of intensive researches, see for example, [START_REF] Lasserre | Nonlinear optimal control: Numerical approximations via moments and lmi-relaxations[END_REF], [START_REF] Chernousko | Control of nonlinear dynamical systems: methods and applications[END_REF], [START_REF] Dinuzzo | Higher order sliding mode controllers with optimal reaching[END_REF], [START_REF] Boltyanski | The Robust Maximum Principle: Theory and Applications[END_REF].

The "bang-bang" control is usual solution of the minimum time control problem. Sometimes, continuous control strategies are preferable in practice, see, for example, [START_REF] Bernstein | Optimal nonlinear, but continuous, feedback control of systems with saturating actuators[END_REF], [START_REF] Glizer | Robust Controllability of Linear Systems[END_REF]. Suboptimal robust control solutions have a great practical interest [START_REF] Bertsekas | Dynamic programming and suboptimal control: A survey from adp to mpc[END_REF], [START_REF] Bohl | Linear feedback vs. time optimal control[END_REF]), since they are usually obtained under additional restrictions to the control forms, which are required for a practical control application and robustness of the control strategy.

The present paper considers the problem of time optimal control synthesis for a chain of integrators controlled by linear dynamic feedbacks. The design concept is based on the Implicit Lyapunov Function (ILF) method [START_REF] Korobov | A general approach to synthesis problem[END_REF], [START_REF] Adamy | Soft variable-structure controls: a survey[END_REF], [START_REF] Polyakov | Finite-time stabilization using implicit lyapunov function technique[END_REF]) in order to present the suboptimal solution for minimum time control problem. This method uses Lyapunov function defined in the implicit form by means of an algebraic equation. Convergence and stability analysis in this case does not require solving of this equation, since the Implicit function theorem (see, for example, [START_REF] Courant | Introduction to calculus and analysis[END_REF]) helps to provide all required analysis using the algebraic equation only. The ILF control design scheme has representation in the form of Linear Matrix Inequalities (LMI), see [START_REF] Polyakov | Finite-time stabilization using implicit lyapunov function technique[END_REF]. The time optimization problem in this case becomes SDP (semidefinite programming) formalization.

The paper is organized as follows. The next section introduces some notations used in the paper. The section 3 presents the problem statement. The section 4 gives some preliminary remarks about finite-time ILF method. Then the control design schemes are given and studied. The aspects of the practical implementation of the developed control schemes are discussed in the section 6. Finally, numerical examples and conclusions are presented.

NOTATION

Through the paper the following notations will be used:

• R + = {x ∈ R : x > 0}, R -= {x ∈ R : x < 0}, where
R is the set of real number; • i = √ -1 is the imaginary unit; • for a differential equation numbered as (.), the time derivative of a function V along the solution of (.) is denoted by dV dt (.) ;

• • is the Euclidian norm in R n , i.e. x = x 2 1 + . . . + x 2 n for x = (x 1 , . . . , x n ) T ∈ R n ; • diag{λ i } n
i=1 is the diagonal matrix with the elements λ i on the main diagonal; • a continuous function σ : R + → R + belongs to the class K if it is monotone increasing and σ(0) = 0; • for a matrix P ∈ R n×n , which has the real spectrum, the minimal and maximal eigenvalues are denoted by λ min (P ) and λ max (P ), respectively; • if P ∈ R n×n then the inequality P > 0 (P ≥ 0, P < 0, P ≤ 0) means that P is symmetric and positive definite (positive semidefinite, negative definite, negative semidefinite).

• C k [a,b]

PROBLEM STATEMENT

Let us consider a single input control system of the form

ẋ = Ax + bu, t ∈ [0, T ], x(0) = x 0 ∈ R n , (1) where T ∈ R + is a finite instant of time, x ∈ R n is the state vector, u ∈ R is a control input, A =      0 1 0 ... 0 0 0 1 ... 0 ... ... ... ... ... 0 0 0 ... 1 0 0 0 ... 0      and b =      0 0 ... 0 1      .
The system (1) describes a controlled chain of integrators. The model of the control system like (1) is motivated by many mechanical and electromechanical applications, see, for example, [START_REF] Chernousko | Control of nonlinear dynamical systems: methods and applications[END_REF], [START_REF] Utkin | Sliding Mode Control in Electro-Mechanical Systems[END_REF], [START_REF] Biagiotti | Time-optimal regulation of a chain of integrators with saturated input and internal variables: An application to trajectory planning[END_REF].

Let us consider the time-optimal control problem

T → min (2) subject to u(•) ∈ {u(•) ∈ L [0,T ] : |u(t)| ≤ u 0 , t ∈ [0, T ]}, x(•) ∈ C 1 [0,T ] : ẋ(t) = Ax(t) + bu(t), x(0) = x 0 , x(T ) = 0, (3)
where x 0 ∈ R n and u 0 ∈ R + are given.

The optimization problem (2), ( 3) is the classical timeoptimal control problem. According to the celebrated theorem by [START_REF] Feldbaum | Optimal processes in systems of automatic control[END_REF] the optimal solution of this problem is the so-called bang-bang control, which is a piecewise constant function u opt (t) ∈ {-u 0 , u 0 } with n switching instants (including t = 0). This paper addresses the suboptimal control design to the optimization problem (2), (3) under some additional assumptions. The first one is the following:

u(•) ∈ U, U = u(•) ∈ C [0,T ] ∩ C ∞ [0,T ) : |u(t)| ≤ u 0 , t ∈ [0, T ] (4)
where u 0 ∈ R + is some given number. In other words, the scalar control input u ∈ R is assumed to be bounded and smooth.

The saturation and smoothness of the input signal are the natural practical restrictions of control systems considered in many papers, see for example, Wang et al.

[2012], [START_REF] Tarbouriech | Control of uncertain systems with bounded inputs[END_REF]. Frequently, these conditions are necessary for a physical realization of a control algorithm.

Since the linear feedback strategies may provide a nearly time-optimal solutions [START_REF] Bohl | Linear feedback vs. time optimal control[END_REF]), then let us restrict additionally admissible control strategies to the class of linear time-varying feedbacks.

The main aim of this paper is to propose an algorithm for time suboptimal control design

u(•) = w T (•)x(•) ∈ U, w = (w 1 , w 2 , ..., w n ) T ∈ R n , w i (•) ∈ C ∞ [0,T ) , (5)
which admits the simple scheme of parameters tuning based on LMI formalism. Such approach allows us to reduce (relax) the infinite dimensional optimization problem (2), (3) to a finite dimensional one, which provides a suboptimal solution.

The last additional assumption is not constructive from the mathematical point of view. However, it is very impor-tant for control practice, which frequently looks for simple solutions for improvement of the control performance.

Another problem, which is discussed in the paper, is the robustness of the proposed control scheme. Robustification of the optimal solutions is a very important task (Boltyanski and Poznyak [2012], [START_REF] Chernousko | Control of nonlinear dynamical systems: methods and applications[END_REF], [START_REF] Dinuzzo | Higher order sliding mode controllers with optimal reaching[END_REF]).

FINITE-TIME ANALYSIS USING IMPLICIT LYAPUNOV FUNCTION

Consider the system of the form ẋ = f (t, x), x(0) = x 0 , (6) where x ∈ R n is the state vector, f : R + × R n → R n is a continuous nonlinear vector field, which is locally Lipschitz outside the origin.

Assume that the origin is an equilibrium point of the system (6), i.e. f (t, 0) = 0. Definition 1. [START_REF] Roxin | On finite stability in control systems[END_REF], [START_REF] Bhat | Finite-time stability of continuous autonomous systems[END_REF]). The origin of the system ( 6) is said to be finite-time stable if for a set V ⊂ R n :

(1) Finite-time attractivity: there exists a function T :

V \ {0} → R + , such that for all x 0 ∈ V \ {0}, x(t, x 0 ) is defined on [0, T (x 0 )) and lim t→T (x0)

x(t, x 0 ) = 0.

(2) Lyapunov stability: there exists a function δ ∈ K such that for all x 0 ∈ V, x(t, x 0 ) ≤ δ( x 0 ).

The function T (•) from Definition 1 is called the settlingtime function of the system (6).

If the set V in Definition 1 coincides with R n then the origin of the system (6) is globally finite time stable.

The next theorem presents the extension of the ILF method [START_REF] Adamy | Soft variable-structure controls: a survey[END_REF] for finite-time stability analysis. Theorem 2. [START_REF] Polyakov | Finite-time stabilization using implicit lyapunov function technique[END_REF]). If there exists a continuous function

Q : R + × R n → R (V, x) → Q(V, x) such that C1 ) Q is continuously differentiable ∀x ∈ R n \{0} and ∀V ∈ R + ; C2 ) for any x ∈ R n \{0} there exist V -∈ R + and V + ∈ R + : Q(V -, x) < 0 < Q(V + , x); C3 ) lim x→0 (V,x)∈Ω V = 0 + , lim V →0 + (V,x)∈Ω x = 0, lim x →∞ (V,x)∈Ω V = +∞, where Ω = (V, x) ∈ R n+1 : Q(V, x) = 0 ; C4 ) the inequality -∞ < ∂Q(V,x) ∂V < 0 holds ∀V ∈ R + and ∀x ∈ R n \{0}; C5 ) for ∀t ∈ R + and ∀(V, x) ∈ Ω; ∂Q(V, x) ∂x f (t, x) ≤ cV 1-µ ∂Q(V, x)
∂V where c > 0 and 0 < µ ≤ 1 are some constants, then the origin of system ( 6) is globally finite-time stable with the settling-time estimate T (x 0 ) ≤

V µ 0 cµ , where V 0 ∈ R + : Q(V 0 , x 0 ) = 0.
The ILF method defines a Lyapunov function V for the a system in the implicit form, for instance, as a solution of the algebraic equation Q(V, x) = 0. Conditions C1 )-C4 ) of Theorem 2 guarantee the properness (positive definiteness, radial unboundedness) of the corresponding solution. Due to the classical implicit function theorem (see, for example, [START_REF] Courant | Introduction to calculus and analysis[END_REF]) the partial derivative of the Lyapunov function V can be calculated as

follows ∂V ∂x = -∂Q ∂V -1 ∂Q ∂x .
Therefore, the condition C5 ) guarantees negative definiteness of the total derivative of V and finite-time stability of the system (6).

CONTROL DESIGN

Finite-Time Control

Introduce the Implicit Lyapunov Function

Q(V, x) := x ⊤ D µ (V -1 )P D µ (V -1 )x -1, (7)
where D µ (λ) is the dilation matrix of the form

D µ (λ) = diag{λ 1+(n-i)µ } n i=1
, 0 < µ ≤ 1 and P ∈ R n×n is a symmetric positive definite matrix, i.e. P = P ⊤ > 0.

The following Theorem presents the main scheme of the ILF control design. It refines the result of [START_REF] Polyakov | Finite-time stabilization using implicit lyapunov function technique[END_REF] providing the precise estimate of the settling time. Denote H µ := diag{1 + (n -i)µ} n i=1 . Theorem 3. Let X ∈ R n×n and y ∈ R 1×n satisfy the following system of LMIs:

AX + XA ⊤ + by + y ⊤ b ⊤ + H µ X + XH µ = 0, XH µ + H µ X > 0, X > 0, (8)
for some µ ∈ (0, 1]. Then the control of the form

u(V, x) = V 1-µ kD µ (V -1 )x, (9)
where k = yX -1 , V ∈ R + : Q(V, x) = 0 and Q(V, x) is defined by ( 7) with P = X -1 , stabilizes the system (1) in a finite time with settling-time function defined as follows

T (x 0 ) = V µ 0 µ , (10)
where

V 0 ∈ R + : Q(V 0 , x 0 ) = 0.
Proof. The function Q(V, x) defined by ( 7) satisfies the conditions C1 )-C3 ) of Theorem 2. Indeed, it is continuously differentiable for all V ∈ R + and ∀x ∈ R n . Since P > 0 then the following chain of inequalities λ min (P )

x 2 max{V 1+(n-1)µ , V } ≤ Q(V, x) + 1 ≤ λ max (P ) x 2 min{V 1+(n-1)µ , V } implies that for any x ∈ R n \{0} there exist V -∈ R + and V + ∈ R + : Q(V -, x) < 0 < Q(V + , x). Moreover, if Q(V, x) = 0 then the same chain of inequalities gives min{V 1+(n-1)µ , V } λ max (P ) ≤ x 2 ≤ max{V 1+(n-1)µ , V } λ min (P ) .
It follows that the condition C3 ) of Theorem 2 holds.

Since

∂Q ∂V = -V -1 x ⊤ D µ (V -1 )(H µ P + P H µ )D µ (V -1 )x
then (8) and P := X -1 implies H µ P + P H µ > 0 and ∂Q ∂V < 0 for ∀V ∈ R + and x ∈ R n \{0}. So the condition C4 ) of Theorem 2 also holds. In this case we have ∂Q ∂x (Ax + bu(x)) = 2x ⊤ D r (V -1 )P D µ (V -1 )(Ax + bu(x)).

Taking into account that D r (V -1 )AD

-1 r (V -1 ) = V -µ A and D µ (V -1 )bu(x) = V -µ bkD µ (V -1 )x we obtain ∂Q ∂x (Ax + bu(x)) = V -µ x ⊤ D µ (V -1 ) P (A + bk) + (A + bk) ⊤ P D µ (V -1 )x Therefore dV dt (1) = - ∂Q ∂V -1 ∂Q ∂x (Ax + bu(x)) = x ⊤ D µ (V -1 ) P (A + bk) + (A + bk) ⊤ P D µ (V -1 )x x ⊤ D µ (V -1 )(H µ P + P H µ )D µ (V -1 )x V 1-µ = z ⊤ P (A + bk) + (A + bk) ⊤ P z z ⊤ (H µ P + P H µ )z V 1-µ = -V 1-µ
where z := z(V, x) = D µ (V -1 )x. The last equality implies the finite-time stability of the system (1) and the representation (10) for the settling-time function.

The main advantage of the presented control scheme is the LMI representation of the procedure for parameters tuning. Any existing LMI solver can be used in order to select the control parameters. Proposition 4. [START_REF] Polyakov | Robust nonasymptotic stabilization of linear systems[END_REF]). The LMI system (8) is feasible for all µ > 0.

A practical implementation of the control (9) admits a simple realization if the initial state x(0) = x 0 ∈ R n of the system (1) is given. Indeed, since the function Q(V 0 , x 0 ) satisfies the conditions of Theorem 2 then the equation Q(V 0 , x 0 ) = 0 has a unique positive solution V 0 ∈ R + for any given x 0 ∈ R n . The corresponding solution can be found numerically using, for example, the simplest bisection method.

Since the function V is a solution of (see the proof of Theorem 3) the Cauchy problem V (t) = -V 1-µ , V (0) = V 0 , then it can be found analytically

V (t) = (V µ 0 -µt) 1/µ , t ∈ [0, V µ 0 /µ].
(11) If the linear dynamic feedback u(•) has the form (9) with V = V (t), where V (t) is defined by ( 11) with V 0 satisfying the equation Q(V 0 , x 0 ) = 0 for a given x 0 ∈ R n , then the condition x(V µ 0 /µ) = 0 fulfills on the corresponding trajectory of the closed system (1).

Suboptimal Finite-Time Control

The value T = V µ 0 /µ has to be minimized under the constraint u(•) ∈ U, where u is defined by (9), in order to obtain a suboptimal solution to the problem (2), (3). Let us formulate a finite dimensional optimization problem with the LMI constraints, which provides the required suboptimal solution. Corollary 5. If for some µ ∈ (0, 1], x 0 ∈ R n \{0} and T ∈ R + the LMI system (8) is feasible together with 1

x ⊤ 0 x 0 D µ (µT ) 1/µ XD µ (µT ) 1/µ ≥ 0, (12)

  is a set of functions defined on [a, b], which are continuously differentiable at least up to the order k. • L [a,b] is a set of Lebesgue integrable functions defined on [a, b].

where X ∈ R n×n and y ∈ R 1×n , then the ILF control (9) guarantees that the settling time function T (x 0 ) of the system (1) satisfies the inequality T (x 0 ) ≤ T , where x 0 ∈ R n is a given initial state.

Proof. Denote V 0 ∈ R + : Q(V 0 , x 0 ) = 0.

The LMI (12) implies

According to Theorem 3 we obtain T (x 0 ) =

Minimization of the parameter T ∈ R + under the LMI constraints ( 8), ( 12) minimizes the settling time of the closed-loop system (1),( 9). If the LMI ( 8) is feasible then the settling time T (x 0 ) can be made infinitely small if the magnitude of the control signal is not bounded. Indeed, let the (y 0 , X 0 ) satisfies the LMI (8), then y = αy 0 , X = αX 0 is also satisfies (8) for any α ∈ R + . Since X 0 > 0 then for any T > 0 and any x 0 ∈ R n \{0} there exists α ∈ R + such that the LMI (12) holds. Corollary 6. Let for some µ ∈ (0, 1], T ∈ R + and u 0 ∈ R + the LMI system (8), ( 12) is feasible together with

where X ∈ R n×n and y ∈ R 1×n . Then the ILF control ( 9) is bounded by u 0 along the generated trajectory x(t) of the closed-loop system, i.e. |u(V (t),

Proof. Using Shur complement the LMI (13) can be rewritten in the form (µT ) 2-2µ µ u 2 0 y ⊤ y ≤ X. Taking into account X = P -1 and y = kP -1 we obtain (µT )

]. The proven corollary allows designing the ILF control of the given maximum magnitude. If µ = 1 then the inequality (13) becomes independent of T and the control is globally bounded. It is important to remark that the design procedure is again formulated in the LMI form.

For any fixed µ and T the matrix inequalities ( 8), ( 12), ( 13) become LMIs, which can be easily solved using any LMI solver (for example, SeDuMi).

Summarizing the obtained results we can formulate the following finite dimensional optimization problem for a fixed µ ∈ (0, 1]:

T → min

X∈R n ,y∈R 1×n

subject to ( 8), ( 12), ( 13). ( 14)

For the fixed parameter µ the presented optimization problem can be easily solved using, for example, the method of intervals (see, for example, [START_REF] Tuan | Advanced global optimization algorithms for parameterized lmis[END_REF]).

Minimization with respect to µ can be done applying any derivative free method (for example, the function fminsearch of MATLAB).

ON ROBUST DIGITAL IMPLEMENTATION OF SUBOPTIMAL ILF CONTROL

As it was remarked above the value of the Lyapunov function can be calculated along the trajectory as follows

The suboptimal linear dynamic feedback u(•) can be defined in the form (9) with V = V (t), where V (t) is defined by ( 11) and by the equation Q(V 0 , x 0 ) = 0 for a given x 0 ∈ R n . In this case the trajectory of the closed system (1) satisfies the condition x(V µ 0 /µ) = 0. In practice the original system (1) has some disturbances and uncertainties that can reduce the quality of the presented dynamic control. Therefore, some practical robustification procedures are required for the suboptimal ILF control scheme.

Sampled-time realization of the ILF control

The ideas of the sampled-time realization of the ILF control are discussed in [START_REF] Polyakov | Finite-time stabilization using implicit lyapunov function technique[END_REF]. They are mainly based on two following corollaries, which can also be proven for the ILF control defined in Theorem 3. Corollary 7. [START_REF] Polyakov | Finite-time stabilization using implicit lyapunov function technique[END_REF]) Let the control u(V, x) be defined according to Theorem 3 then the control u 0 (x) = u(V 0 , x) is the linear stabilizing feedback control for the system (1) for any given V 0 ∈ R + .

The next corollary will help us to analyze the discrete-time version of the developed control schemes. Corollary 8. [START_REF] Polyakov | Finite-time stabilization using implicit lyapunov function technique[END_REF]) Let {t i } ∞ i=0 be an arbitrary strictly increasing sequence of time instants, 0 = t 0 < t 1 < t 2 < .... Let the function u(V, x) is defined according to Theorem 3 then the origin of the system (1) with the switching control u(x)

The last corollary provides that the sampled-time realization of the ILF control keeps the stability property of the closed-loop system (1) independently on the sampling period.

Let the control u(V, x) of the form (9) be realized in a digital device and the parameter V may change its value at some time instants: t 0 = 0, t i > 0, i = 1, 2, .... Denote V i := V (t i ) and x i := x(t i ). On the time interval [t i , t i+1 ) the control has the form u(V i , x), which is a linear stabilizing feedback for any V i ∈ R + (see, Corollary 7).

The simplest scheme for selection of the switching control parameter V i is described by the following algorithm.

A L G O R I T H M [START_REF] Polyakov | Finite-time stabilization using implicit lyapunov function technique[END_REF])

Remark, STEP of Algorithm 6.1 can be realized just once in each sampled time instant.

The parameter V min defines lower possible value of V . In practice, this parameter cannot be selected arbitrary small due to finite numerical precision of digital devices.

Adaptive scheme

An adaptive scheme can also be presented for practical implementation of ILF control.

Theorem 3 proves the following dynamic equation for the

If this equation is calculated on-line with the initial condition V (0) = V 0 > 0 : x(0) T D µ (V -1 0 )P D µ (V -1 0 )x(0) = 1 then the corresponding trajectory x(t) of the system (1) with the control ũ(V (t), x(t)) converges to the origin in a finite time (see, Theorem 2 and the previous section).

In order to provide some robustness to the presented control scheme the dynamic equation V = -V 1-µ can be modified as follows

where H[ρ] is a Heaviside step function

The modified dynamic equation ( 16) guarantees that the function V (t) will be decreasing only in the case

if s belongs to the ellipsoidal level set of the ILF. Otherwise, V = 0 and the function V (t) is constant. Since ũ(V, s) is a linear stabilizing feedback for any fixed V > 0 then the adaptive scheme (16) of the ILF implementation will guarantee at least asymptotic stabilization of the system (1). For the chattering reduction the function V (t) can also be restricted from below by some minimal value V min .

The Heaviside function in the equation ( 16) can be replaced with sign[-Q(V (t), x(t)] providing additional computational robustness of the digital ILF control application. A more detailed study of the ILF control implementation goes out of the scope of this paper. It is suggested as a topic for future research. The numerical simulations given below demonstrate an effectiveness of the presented simplest algorithms of ILF control implementation.

NUMERICAL EXAMPLE

Consider the system (1) for n = 2 under control restriction |u| ≤ 1.

It is well known, see, for example, [START_REF] Chernousko | Control of nonlinear dynamical systems: methods and applications[END_REF], that the time optimal feedback law for this case has the discontinuous form

The continuous stabilizing feedback (9) designed following the optimization procedure ( 14) for µ = 1 and x 0 = (1, 0) has the parameters P = 49.6139 18.8965 18.8965 9.4482 , k = (-5.2511 -3.0000).

The numerical solving of ODE for the closed-loop system have been done by the Euler method with the fixed step size and the discrete time application of the finite time ILF control is realized by the scheme presented in Algorithm 6.1. The adaptive scheme has also been verified by numerical simulation. It provides very similar results.

The results of numerical simulations for h = 0.001 and for V min = 0.001 are shown on Fig. 1. It was expectable that the optimal controller provides the faster transition. However, if the control will be applied with larger sampling period, which frequently takes a place in real applications, then the situation is changing. Fig. 2 depicts the simulation results for h = 0.05 and V min = 0.05. Such defect of the optimal controller is known as "chattering phenomenon" (see, [START_REF] Utkin | Sliding Mode Control in Electro-Mechanical Systems[END_REF]), which appears in time-optimal solutions realized by sliding mode control feedbacks [START_REF] Dinuzzo | Higher order sliding mode controllers with optimal reaching[END_REF]).

CONCLUSIONS

In the paper a new suboptimal control algorithm is presented. It has the following advantages:


Julien Toulouse
email: julien.toulouse@upmc.fr

Iann C Gerber

Georg Jansen

Andreas Savin

János G Ángyán
email: janos.angyan@crm2.uhp-nancy.fr

Adiabatic-connection fluctuation-dissipation density-functional theory based on range separation

An adiabatic-connection fluctuation-dissipation theorem approach based on a range separation of electron-electron interactions is proposed. It involves a rigorous combination of short-range density functional and long-range random phase approximations. This method corrects several shortcomings of the standard random phase approximation and it is particularly well suited for describing weaklybound van der Waals systems, as demonstrated on the challenging cases of the dimers Be2 and Ne2.

Density functional theory (DFT) is a powerful approach for electronic-structure calculations of molecular and condensed-matter systems [1]. However, one difficulty in its Kohn-Sham (KS) formulation using local density or generalized-gradient approximations (LDA and GGA) is the description of non-local correlation effects, such as those involved in weak van der Waals complexes, bound by London dispersion forces [2]. The adiabaticconnection fluctuation-dissipation theorem (ACFDT) approach is one of the most promising ways of constructing highly non-local correlation functionals. This approach, introduced in wave function theory [3] and in DFT [4,5], consists in extracting non-local ground-state correlations from the linear charge density response function.

Recently, the ACFDT approach has received renewed interest for implementing the random phase approximation (RPA) or other related approximations for atoms, molecules and solids [6][7][8][9][10][11]. The RPA correlation energy is consistent with the use of the exact, self-interactionfree exchange energy. In spite of a number of encouraging results, such as the correct description of dispersion forces at large separation [12], the proper reproduction of cohesive energies and lattice constants of solids [10,13,14] and an improved description of bond dissociation [6,7,15], several aspects of the RPA are still unsatisfactory.

First, the RPA is a poor approximation to short-range correlations, leading to correlation energies that are far too negative [16]. Second, in a Gaussian localized basis, RPA calculations have a slow convergence with respect to the basis size [6]. Third, the presence of an unphysical maximum (bump) at medium distances in dissociation curves of simple diatomic molecules [6,15] indicates an inherent problem which has not yet a fully clarified origin. Fourth, although in principle the orbitals should be calculated self-consistently [17], most RPA implementations consist of a post-KS single-iteration calculation, making the choice of the input orbitals sometimes critical. Last but not least, although the main advantage of the RPA is supposed to be the description of dispersion forces, rare gas dimer potential curves calculated from LDA or GGA orbitals are often qualitatively wrong, as shown later.

The poor short-range behavior can be corrected by adding a GGA functional constructed from the difference of the exact and RPA correlation energies of the uniform electron gas [16], but this so-called RPA+ technique does not lead to consistent improvement [6]. One can go beyond RPA by including exchange-correlation (xc) kernels [8,9,18], but so far it remains imperfect, e.g. local xc kernels produce pair densities that diverge at small interparticle distances [9,18].

In a similar spirit as in the work of Kohn et al. [19], we propose an ACFDT approach based on a range separation of electron-electron interactions. It involves a rigorous combination of a short-range density functional with one of the possible long-range generalizations of the RPA. The method offers a solution for several of the aforementioned difficulties and is particularly well suited for the description of weakly-bound van der Waals systems.

Theory. In the range-separated multideterminant extension of the KS scheme, an alternative approach to DFT (see, e.g., Ref. 20), the exact ground-state energy of an electronic system is expressed as

E = min Ψ Ψ| T + Vne + Ŵ lr ee |Ψ + E sr Hxc [n Ψ ] , (1)
where T is the kinetic energy operator, Vne is the nuclei-electron interaction operator, Ŵ lr ee = (1/2) dr 1 dr 2 w lr ee (r 12 )n 2 (r 1 , r 2 ) is a longrange electron-electron interaction written with w lr ee (r) = erf(µr)/r and the pair-density operator n2 (r 1 , r 2 ) = n(r 1 )n(r 2 ) -n(r 1 )δ(r 1r 2 ), and E sr Hxc [n] is the corresponding µ-dependent short-range Hartreeexchange-correlation (Hxc) density functional that Eq. ( 1) defines. The minimizing multideterminant wave function, denoted by Ψ lr , corresponds to a long-range interacting effective Hamiltonian and yields the exact density. The parameter µ in the error function controls the range of the separation. For µ = 0, the standard KS scheme is recovered: w lr ee (r) vanishes, Ψ lr reduces to the non-interacting KS wave function and E sr Hxc [n] becomes the usual Hxc functional. For µ → ∞, the usual wave function formulation of the electronic problem is retrieved: w lr ee (r) becomes the full Coulomb interaction, E sr Hxc [n] vanishes and Ψ lr becomes the exact groundstate wave function. For intermediate values of µ, the interaction effects are divided between the long-range interacting wave function Ψ lr and the short-range density functional E sr Hxc [n], and one expects to find better approximations for each piece. Short-range LDA [21,22] and several beyond-LDA approximations [20,[23][24][25] have been proposed for E sr Hxc . Here, we use an ACFDT approach for the long-range part of the calculation.

In a first step, the minimization in Eq. ( 1) is restricted to single-determinant wave functions Φ, leading to a range-separated hybrid (RSH) scheme [26]

E RSH = min Φ Φ| T + Vne + Ŵ lr ee |Φ + E sr Hxc [n Φ ] , (2)
which, in contrast to some range-separated KS schemes [27,28], does not include long-range correlation. The minimizing determinant Φ 0 is given by the self-consistent Euler-Lagrange equation

T + Vne + V lr Hx,HF + V sr Hxc |Φ 0 = E 0 |Φ 0 , (3)
where V lr Hx,HF is a Hartree-Fock (HF) type longrange Hartree-exchange (Hx) potential, V sr Hxc drδE sr Hxc [n Φ0 ]/δn(r)n(r) is the short-range local Hxc potential and E 0 is the Lagrange multiplier for the normalization constraint. As usual, V lr Hx,HF is the sum of a local Hartree part V lr,µ H and a non-local exchange part V lr x,HF . The RSH scheme does not yield the exact energy and density, even with the exact short-range functional E sr Hxc . Nevertheless, the RSH approximation can be used as a reference to express the exact energy as

E = E RSH + E lr c , (4)
defining the long-range correlation energy E lr c , for which we will now give an adiabatic connection formula. For that, we introduce the following energy expression with a formal coupling constant λ

E λ = min Ψ Ψ| T + Vne + V lr Hx,HF + λ Ŵ lr |Ψ +E sr Hxc [n Ψ ] , (5)
where Ψ is a multideterminant wave function, Ŵ lr is the long-range fluctuation potential operator

Ŵ lr = Ŵ lr ee -V lr Hx,HF , (6)
and E sr Hxc is the previously-defined λ-independent shortrange Hxc functional. The minimizing wave function is denoted by Ψ lr λ . For λ = 1, the physical energy E = E λ=1 and density are recovered, as Eq. ( 5) reduces to Eq. (1).

For λ = 0, the minimizing wave function is the RSH determinant Ψ lr λ=0 = Φ 0 . Note that, because the density at λ = 0 is not exact, the density is supposed to vary along this adiabatic connection. Taking the derivative of E λ with respect to λ, noting that E λ is stationary with respect to Ψ lr λ , and reintegrating between λ = 0 and λ = 1 gives

E = E λ=0 + 1 0 dλ Ψ lr λ | Ŵ lr |Ψ lr λ , (7)
with

E λ=0 = Φ 0 | T+ Vne + V lr Hx,HF |Φ 0 +E sr Hxc [n Φ0 ] = E RSH - Φ 0 | Ŵ lr |Φ 0 .
Thus, the long-range correlation energy is

E lr c = 1 0 dλ Ψ lr λ | Ŵ lr |Ψ lr λ -Φ 0 | Ŵ lr |Φ 0 , (8)
or, using a compact notation,

E lr c = 1 2 1 0 dλ Tr w lr * P lr c,λ , (9)
where w lr and P lr c,λ are four-index representations of the fluctuation potential and the correlation contribution of the two-particle density matrix in a one-electron basis, * stands for contraction of two indices and Tr is the trace over the remaining two indices. The fluctuationdissipation theorem is then used to express P lr c,λ with the imaginary-frequency four-point polarizability χ lr λ (iu) corresponding to the wave function Ψ lr λ (see, e.g., Ref.

P lr c,λ = - 1 2π ∞ -∞ du e -u0 + χ lr λ (iu) -χ 0 (iu) +∆ lr λ , 3)
where χ 0 (iu) is the four-point polarizability for the RSH effective Hamiltonian of Eq. ( 3) and ∆ lr λ is the contribution coming from the variation of the one-particle density-matrix along the adiabatic connection. The expression of ∆ lr λ is straightforward but it is sufficient to write it as ∆ lr

λ = F [G lr λ ] -F [G 0 ]
where F is a known functional, G lr λ is the two-point one-particle Green function corresponding to the wave function Ψ lr λ and G 0 is the two-point RSH Green function. Along the adiabatic connection of Eq. ( 5), the Green function G lr λ satisfies a self-consistent Dyson equation

G lr λ -1 = G -1 0 -λ Σ lr Hx [G lr λ ] -Σ lr Hx [G 0 ] -Σ lr c,λ [G lr λ ], (11)
where λΣ lr Hx and Σ lr c,λ are the Hartree-exchange and correlation self-energies associated with the long-range interaction w lr ee . The long-range polarizability is given by the solution of the Bethe-Salpeter-type equation (see Ref. 29)

χ lr λ -1 = χ lr IP,λ -1 -λf lr Hx -f lr c,λ , (12)
where χ lr IP,λ is an independent-particle (IP) polarizability whose expression is a frequency convolution of two Green functions G lr λ , abbreviated as χ lr IP,λ = -iG lr λ G lr λ , and λf lr Hx = iλδΣ lr Hx /δG lr λ and f lr c,λ = iδΣ lr c,λ /δG lr λ are longrange HF-type and correlation kernels.

So far, the theory is in principle exact. In the following we introduce the approximation

Σ lr c,λ = 0, ( 13
)
which corresponds to neglecting long-range correlation in all one-electron properties. Indeed, from Eq. ( 11), this approximation implies that the Green function remains unchanged along the adiabatic connection, i.e. G lr λ = G 0 and thus ∆ lr λ = 0. It also follows that f lr c,λ = 0 and χ lr IP,λ = χ 0 , and Eq. ( 12) then has the structure of the RPA with HF exchange kernel (sometimes refers to as linear response time-dependent Hartree-Fock theory or full RPA) that we will designate here by RPAx, by opposition to standard RPA, without exchange kernel (sometimes also called direct RPA). In the basis of RSH spatial orbitals, for spin-restricted closed-shell systems, the long-range correlation energy then becomes

E lr c = 1 2 1 0 dλ iajb ij| ŵlr ee |ab P lr c,λ iajb , (14)
where ia and jb refer to excitations from occupied (i or j) to virtual (a or b) orbitals, ij| ŵlr ee |ab are the two-electron integrals with long-range interaction, and (P lr c,λ ) iajb are the matrix elements of the spin-singletadapted P lr c,λ . The one-electron term V lr Hx,HF in Eq. ( 6) does not contribute to E lr c because of the occupied-virtual structure of P lr c,λ . Only singlet excitations contribute to Eq. ( 14), since the two-electron integrals involved vanish for triplet excitations. Note that alternative (but inequivalent) RPAx correlation energy expressions, such as the plasmon formula of Ref. 3 and the closely related ring CCD approximation of Ref. 30, require contributions from both singlet and triplet excitations, which may be problematic in systems displaying triplet instabilities, such as Be 2 . Following the technique proposed by Furche [6], P lr c,λ can be obtained as

P lr c,λ = 2 (A λ -B λ ) 1/2 M -1/2 λ (A λ -B λ ) 1/2 -1 ,( 15
)
with

M λ = (A λ -B λ ) 1/2 (A λ + B λ )(A λ -B λ )
where ǫ i are the RSH orbital eigenvalues. This method will be referred to as RSH+RPAx. In the limit of µ = 0, it reduces to the standard KS scheme, while for µ → ∞ it becomes a full-range ACFDT RPAx approach (with HF orbitals). We note that at secondorder in the interaction w lr ee the RSH+RPAx reduces to the RSH+MP2 method [26,31].

Computational details. Equations ( 14)-( 17) have been implemented in the time-dependent DFT development module [32] of MOLPRO 2008 [START_REF] Werner | Molpro, version 2008.2, a package of ab initio programs[END_REF]. We perform a selfconsistent RSH calculation with the short-range PBE xc functional of Ref. 24 and add the long-range RPAx correlation energy calculated with RSH orbitals. The range separation parameter is taken at µ = 0.5, in agreement with previous studies [28], without trying to fit it. The λ-integration in Eq. ( 14) is done by a 7-point Gauss-Legendre quadrature [6]. We use large Dunning basis sets [START_REF] Dunning | [END_REF] and remove the basis set superposition error (BSSE) by the counterpoise method. The full-range RPA and RPAx calculations have been done with PBE [35] and HF orbitals, respectively. The computational cost of RSH+RPAx is essentially identical to that of full-range ACFDT RPA.

Results and discussion. Figure 1 shows RPA, RPAx and RSH+RPAx total energies with respect to the basis size for Ne 2 . In contrast with full-range RPA and RPAx, a fast convergence is observed for RSH+RPAx, similar to that of standard KS calculations. This improvement is explained by the fact that short-range correlations are compactly described by the short-range density functional. The reduced basis dependence of RSH+RPAx also means smaller BSSE [START_REF]For Ne2 at equilibrium distance with aug-cc-pVQZ basis[END_REF]. Another important point illustrated by Fig. 1 is that the large RPA overestimation of the total energy, apparent for a large enough basis, is remedied by the RSH+RPAx method thanks to a more accurate description of short-range correlations.

Figure 2 shows the interaction energy curves of Be 2 and Ne 2 . The RPA (with PBE orbitals) fails badly: a large bump for Be 2 and an almost completely repulsive curve for Ne 2 are observed. A spectacular improvement is obtained with the RSH+RPAx method, which gives physically correct curves, especially accurate at medium and large distances. It also improves over both MP2 and RPAx. It has been verified that the short-range LDA xc functional of Ref. 21 provides very similar RSH+RPAx interaction curves, indicating a low sensitivity of the method with respect to the short-range functional.

The proposed RSH+RPAx method overcomes many limitations of the RPA or related approaches. The results show that it has the potential to describe successfully weakly-bound van der Waals systems at all distances. It is expected to supersede the RSH+MP2 method [26], especially for systems with small electronic gap.

2 FIG. 1 :

 21 FIG.1:(Color online) RPA, RPAx and RSH+RPAx total energies of Ne2 at equilibrium distance (5.84 Bohr) with respect to the basis size.

2 FIG. 2 :

 22 FIG. 2: (Color online) Interaction energy curves of Be2 and Ne2 calculated in RPA, MP2, RPAx and RSH+RPAx. The basis is cc-pV5Z for Be2 and aug-cc-pV5Z for Ne2. The accurate curves are from Ref. 36 and Ref. 37.

  1/2 and the singlet orbital rotation Hessians (A λ ) iajb = (ǫ a -ǫ i )δ ij δ ab +2λ aj| ŵlr ee |ib -λ aj| ŵlr ee |bi , (16) (B λ ) iajb = 2λ ab| ŵlr ee |ij -λ ab| ŵlr ee |ji ,

Acknowledgements. We thank F. Furche, J. F. Dobson, T. Gould and G. E. Scuseria for discussions. This work was supported by ANR (07-BLAN-0272).


O Isnard

E G Moshopoulou

J Prchal
email: prchal@karlov.mff.cuni.cz

P Javorsky

Specific heat measurements and structural investigation of CeCu 6-x Sn x compounds

Keywords: 61.05.F-Neutron diffraction and scattering, 75.40.-s Critical-point effects, specific heats, short-range order, 75.30.Mb Valence fluctuation, Kondo lattice, and heavy-fermion phenomena crystal structure, Kondo interactions, heavy fermion, specific heat

The evolution of the crystal structure and some magnetic properties of the heavy fermion material CeCu 6-x Sn x (x = 0, 0.25, 0.65, 0.75, 0.85 and 1.0) has been studied by powder neutron diffraction and by specific heat measurements. The substitution of Cu by Sn suppresses the temperature induced orthorhombic to monoclinic transition, known to occur in the pure CeCu 6 phase. No structural phase transition has been observed in these samples as a function of x but the cell volume increases considerably in an anisotropic way. Sn occupies preferentially the special Cu crystallographic site which is next to each of the four Ce atoms in the unit cell. The transition to antiferromagnetic order, characterising the samples with higher x, is sensitive to both x and magnetic field. The results are discussed in the context of the competition between Kondo and RKKY interactions in disordered or not heavy-fermion systems and reveal an interesting interplay between composition, structure and magnetism in CeCu 6-x Sn x .

Introduction

CeCu 6 is a well-known heavy-fermion material characterised by an exceptionally high electronic contribution (γ) to the specific heat of about 1.6 J/mol K 2 for single crystals at T = 0 K (ref. 1) and Fermi liquid behaviour. Since its discovery in 1984 (ref. 1) it has been extensively investigated; especially its physical properties were probed down to 250 μK. The material does not superconduct but below 1.5 K it exhibits temperature independent magnetic correlations (ref.

2) and it orders antiferromagnetically at 2 mK (ref. 3). CeCu 6 is characterised by a rather unusual for a heavy-fermion material, structural phase transition from an orthorhombic crystal structure (space group Pnma) at room temperature (RT) to a monoclinic structure (s. g. P2 1 /c) at around 230 K with a monoclinic angle of β ≈ 91.58° (ref. 4,5,6). Since such transition was also found in the nonmagnetic analogue LaCu 6 at much higher though temperature (around 490 K) (ref. 6), it was concluded that electronic effects involving the 4f electrons are not responsible for the structural phase transition, i.e. the lattice degrees of freedom are not related to the heavy-fermion properties of CeCu 6 .

Interestingly, Au-doped CeCu 6 exhibits different behaviour from its parent compound. Since the Au atoms are larger than the Cu ions, the crystalline lattice expands and therefore substitution of Cu by Au in CeCu 6 is equivalent to negative chemical pressure on the lattice. The system CeCu 6-

x Au x exhibits long-range antiferromagnetic order at x c ≤x≤1, where x c = 0.1, below Néel temperature, T N , ranging from ~ 0 K (for x c = 0.1) to 2.3 K (for x = 1) (ref. 7). Especially the material CeCu 5.9 Au 0.1 has been extensively investigated over the past twenty years because, while it is metallic, it exhibits non-Fermi liquid properties: linear resistivity between 0.2 and 0.8 K, logarithmic specific heat between 0.08 and 1 K and power-law susceptibility between 0.08 and 3 K (ref. 8). Generally (as indeed has been observed for CeCu 6-x Au x ), such non-Fermi liquid regime is related to a quantum critical point, i.e. to a 2 nd order phase transition typically from paramagnetism to antiferromagnetism at around T = 0 K as a function of a control parameter such as composition, pressure or magnetic field. Among the most common factors that are currently considered to influence the appearance of such non-Fermi liquid behaviour are disorder, magnetic anisotropy, strong Kondo screening. Interestingly, in the case of CeCu 6-x Au x , it was found (ref. 9) that the structural distortion (discussed previously for the parent compound CeCu 6 ) competes with antiferromagnetism, i.e. the structural distortion and the long-range antiferromagnetic order disappear at exactly the same critical concentration x c = 0.13 where non-Fermi liquid behaviour was observed. This observation strongly suggests that the lattice distortion and the concomitant structural fluctuations at x c = 0.13, which are a source of disorder, are also a source of the non-Fermi liquid behaviour of the CeCu 5.9 Au 0.1 . Signatures of such non-Fermi liquid behaviour are found at the temperature dependence of magnetic susceptibility, resistivity and specific heat of heavy-fermion materials. Considerably fewer experimental data exist on the temperature dependence of the thermoelectric coefficients (thermopower, Nerst coefficient) of non-Fermi liquid systems. A few available results reveal unexpectedly enhanced or reduced thermoelectric coefficients. Therefore, the relationship between non-Fermi liquid behaviour, disorder (including lattice distortions and concomitant structural fluctuations) and enhanced or not thermoelectric properties is currently an open question in the physics of heavy-electron materials. The line of research presented here has been conceived as a step in this direction.

Recently, substitution of Cu by Sn has been found to significantly reduce the thermopower of CeCu 6 (ref. 10). Moreover, a crossover from Kondo lattice to antiferromagnetic ordering has been reported to take place in CeCu 6-x Sn x (ref. 10). The phases CeCu 5.25 Sn 0.75 and CeCu 5 Sn order antiferromagnetically at T N = 4 K and 10 K respectively. In order to investigate whether the interplay between structural and magnetic degrees of freedom (discussed above for Au-doped CeCu 6 ) is a generic characteristic of doped-CeCu 6 and gain insight on the factors that are important for its low-temperature transport, magnetic and thermoelectric properties, we investigated the phases CeCu 6-x Sn x for x = 0, 0.25, 0.65, 0.75, 0.85 and 1. Following the earlier report on the properties of CeCu 6-x Sn x phases (ref. 10), we present here an investigation combining a structural analysis by powder neutron diffraction with specific heat measurements. We mostly focus our study to the samples with high Sn-content where an antiferromagnetic order has been reported to take place.

Experimental 2a. Sample preparation

The samples were prepared by melting the starting metals (cerium 99.8%, electrolytic copper 99.99% and tin 99.99%) in an electric arc furnace. The nominal composition of the CeCu 6-

x Sn x samples was: x = 0, 0.25, 0.65, 0.75, 0.85 and 1.0. The obtained materials were annealed at 500°C for 720 h in vacuum. The quality of the samples was checked using conventional powder Xray diffraction.

2b Powder neutron diffraction

The powder neutron diffraction experiments were performed on the D1B instrument (CRG) (ref. 11) at the Institut Laue Langevin (ILL), Grenoble (France). During the neutron diffraction measurements, a cylindrical vanadium sample holder of 6 mm inner diameter was used. The diffraction patterns were recorded over an angular range 2θ of 80° using a multidetector with a step of 0.2° between each of the 400 3 He detection cells. The neutron wavelength was λ = 2.52 Å selected by the (002) Bragg reflection of a pyrolitic graphite monochromator with take off angle of 2θ = 44.2°. The measurements were carried out in the temperature range 10 to 300 K with a step of 7 K using the orange ILL cryostat. For the refinement of the data collected on D1B, a small angular region around 2θ = 72°, containing a peak originating from the vanadium tail of the cryostat, was excluded.

The data were analysed by the Rietveld structure refinement program FULLPROF (ref. [START_REF] Rodriguez-Carvajal | Fullprof Suite[END_REF].

The agreement factors used in this article are defined as χ² = R wp / R exp , R exp = [(N-P+C)/Σ i w i y i ²] 1/2 and R wp = {Σ i w i (y i -y ci )²/ Σ i w i y i ²} 1/2 where y i and y ci are the observed and calculated counts at the i th step of the pattern respectively and w i is the weights of the observations. N-P+C is the number of degrees of freedom (N is the number of points in the pattern, P the number of refined parameters and C the number of constraint functions). The Bragg agreement factor is defined as

R B = Σ k |I k -I ck | / Σ k |I k |

2c. Specific heat measurements

The specific heat of the CeCu 6-x Sn x (x = 0, 0.25, 0.65, 0.75, 0.85 and 1.0) samples was measured using the Physical Property Measurement System (PPMS-14, Quantum Design), installed at DCMP, Charles University, Prague (Czech Republic). The relaxation 2-τ fitting method was used for the measurement. The experiment was performed in the temperature range 0.35 to 40 K and in magnetic fields up to 10 T. The mass of the samples was in the range of 1 to 5 mg. The He-3 option (QD) was installed to reach temperatures down to 350 mK. obtained by X-ray diffraction investigations (ref. 15,16). This atomic position is the site with the highest coordination number: 14 atoms. Unlike the other Cu sites which are surrounded by 12 to 13 near neighbour atoms (3 of which are Ce atoms), the site filled by Sn is surrounded by 14 near neighbours: 10 Cu and 4 Ce atoms. Consequently, the occupation of this site by Sn maximizes the number of Ce-Sn bonds. Furthermore, it is worth to recall that both Sn atoms (like Au) are larger than copper; this size difference may favour such preferential substitution scheme. In order to test this assumption, a calculation of the atomic cell volume of the CeCu 6 type structure has been carried out. The Wigner-Seitz (W.S.) atomic cell volume of the Voronoi-Dirichlet polyhedron around an atomic site has been calculated as described elsewhere (ref. [19][20][21][22] We now turn to the temperature dependence of the crystal structure of the compounds CeCu 5.25 Sn 0.75 and CeCu 5 Sn. Our Rietveld refinements of the powder neutron diffraction data, collected in the temperature range from 4 to 300 K, demonstrate the absence of any structural phase transition for the investigated compositions. Similar behaviour was previously observed for the CeCu 6-x Au x materials with x>0.13 (ref. [START_REF] Grube | [END_REF]25). Therefore, the substitution of Cu by Sn suppresses the orthorhombic to monoclinic transition, known to occur in the pure CeCu 6 phase and in the CeCu 6-x Au x phases with x<0. 13 (ref. 24, 25). The temperature dependence of the lattice parameters below RT is plotted in Figure III for x = 1. As expected, the increase of the temperature from 4 to 300 K leads to an expansion of the unit cell volume of a bit more than 5 Å 3 for CeCu 5 Sn. It is worth to note that no anomalous behaviour is observed in this thermal evolution, thus confirming the absence of any structural phase transition within this temperature range. The thermal expansion of the lattice parameters a and c is very similar, i.e. they both expand of about 0.035 Å upon increasing the temperature from 5 K to RT. The increase of the b parameter is significantly smaller and reaches 0.022 Å. It is worth to mention that in spite of the similar thermal expansion of a and c, observed in Figure III, given that c is much larger than a (about 10 Å and 8 Å respectively), the expansion of the cell is larger in a direction relatively to c direction. Indeed, an increase of 0.42%, 0.43% and 0.33% is observed for a, b and c respectively. This observation demonstrates that the temperature induced expansion of the unit cell, involves all the lattice parameters and it is not isotropic.

Results and discussion

3b. Specific heat measurements

The specific heat C p of CeCu 6-x Sn x (x = 0, 0.25, 0.65, 0.75, 0.85, 1) has been measured at zero magnetic field down to temperature of 0.4 K. The specific heat of CeCu 6 , plotted as C p /T versus T 2 (Fig. IV), exhibits a rapid increase of C p /T below ≈ 7 K which is quite typical of many heavy-fermion materials and is in agreement with the previous studies of CeCu 6 (see for example ref. 1). The sample having low Sn content x = 0.25 exhibits an enhancement of the C p /T value (Fig. V) at low temperatures compared to CeCu 6 . Such behaviour can be suppressed by applying a magnetic field. A linear fit of the low-temperature data C p /T versus T 2 data of CeCu 6 and CeCu 5.75 Sn 0.25 and extrapolation of the C p /T value down to 0 K leads to a γ value of about 2.550(73) J/mole K 2 for CeCu 5.75 Sn 0.25 that is clearly enhanced compared to γ ≈ 1.386(7) J/mole K 2 obtained here for the reference polycrystalline compound CeCu 6 and to γ ≈ 1.6 J/mole K 2 obtained from the literature for the reference single crystalline compound CeCu 6 .

Figure VI shows the zero-field evolution of C p as a function of T in the temperature range 0.4 -20 K of the series CeCu 6-x Sn x (x = 0, 0.65, 0.75, 0.85, 1). For higher Sn concentration x = 0.65, 0.75, 0.85 and 1 an anomaly of C p /T appears, revealing the occurrence of a magnetic order in these compounds, which has been already observed by magnetic susceptibility measurements at earlier studies (ref. 10). We now focus on the broad peak of the C p /T versus T curve observed at T N . This peak is more broad for x = 0.65 and less broad for 0.85. Samples with x = 0.65, 0.75, 0.85 and 1 present magnetic order at T N = 3.5, 3.7, 4.0 and 4.2 K, respectively. This reveals an increase of the antiferromagnetic ordering temperature T N with the increase of the Sn content. Long-range antiferromagnetic order has also been reported to occur in the CeCu 6-x Au x compounds for 0.1≤x≤1 as mentioned in the introduction.

In order to further understand the nature of the strong electron-electron interactions in doped CeCu 6 in general and in CeCu 6-x Sn x in particular, we investigated the influence of an applied magnetic field on the specific heat of CeCu 6-x Sn x samples with high Sn content. The results are plotted in Figures VII, VIII and IX for x = 0.65; 0.75; and 1 respectively. There is a clear decrease of the specific heat at T N caused by the applied magnetic field. This effect is even more pronounced for CeCu 6-x Sn x with x = 0.65 since at 10 T the broad peak at T N is almost suppressed, whereas the sensitivity to the applied magnetic field of the compounds CeCu 6-x Sn x with x = 0.75 and 1 is much weaker. These results reveal that the higher the Sn content in CeCu 6-x Sn x , the weaker the sensitivity of C p to an applied magnetic field. In other words, in an oversimplified picture where disorder is considered to be the same for all samples, the higher the Sn content, the higher has to be the applied magnetic field in order to suppress the magnetic transition.

The evolution of magnetism of CeCu 6-x Sn x as a function of x and/or the applied magnetic field is a result of a delicate balance between various competing interactions which involve its very basic physics, energy scales and crystal chemistry. In an ordered material such as CeCu 6 or CeCu 5 Sn, the exchange interaction J (which couples the local spins to the conduction electrons) and the density of conduction electrons at the Fermi level influence both the Kondo temperature, T K , below which compensation of the local moments by the cloud of conduction electrons occurs, and the Ruderman-Kittel-Kasuya-Yosida (RKKY) temperature, TRKKY, which characterizes the moment-moment interaction in the material. Consequently, in CeCu 6 and in CeCu 5 Sn the exchange interaction J and the density of conduction electrons control the obvious competition between the single-site Kondo screening and inter-site RKKY spin coupling, i.e. between the low-temperature compensated state and the magnetically ordered state of the uncompensated moments. This competition is quite complex because the two energy scales T K and T RKKY have different dependences not only on the density of conduction electrons at the Fermi level but also on ⎢ J ⎢: T K depends exponentially on ⎢ J ⎢ and T RKKY changes as ⎢ J ⎢ 2 . We note that T K = (1/N(E F ))exp(-1/N(E F )J) and T RKKY = -J 2 N(E F ), where N(E F ) is the density of the electrons at the Fermi level.

Turning now to CeCu 6-x Sn x , it is instructive to recall (cf. Section 3a) that the substitution of Cu by Sn in CeCu 6 takes place at a special Cu site which is next to each of the four Ce atoms in the unit cell, maximizing in this way the number of Ce-Sn bonds. Such substitution and the resulting structural changes modify the exchange parameter ⎢ J ⎢, change the density of conduction electrons and introduce disorder. In the CeCu 6-x Sn x phases where disorder effects are quite pronounced, the competition between Kondo screening and RKKY interaction is very complex. It has been demonstrated by several theoretical and experimental studies on a wide variety of heavy-fermion systems that modest disorder leads to modest variations of J which, because it is amplified exponentially, can lead to wide distributions of T K and thus to dramatic changes of the lowtemperature properties of the heavy-fermion material. We believe that, given that the electronic state of Ce does not vary much because of the substitution of Cu by Sn (ref. 10), the evolution of magnetism as a function of x in CeCu 6-x Sn x and its sensitivity or not to the applied magnetic field results from the complex competition between Kondo and RKKY interactions. This competition is influenced by the changes of ⎢ J ⎢ (which in turn are caused by the variations of the interatomic distances, the type of first neighbours of Ce and the disorder) and by the variations of the density of the conduction electrons in the alloyed samples.

Concluding remarks

CeCu 6-x Sn x is an interesting system to investigate the competition between the two basic types of interactions in heavy-fermion systems, the Kondo and RKKY interaction. Our study, focused mostly on the highly doped samples (where antiferromagnetic order takes place), revealed that no structural phase transition occurs in these samples as a function of x or temperature.

However, the substitution of Cu by Sn increases considerably the cell volume (because the Sn atom is bigger than Cu), provides higher concentration of conduction electrons (because Sn has more valence electrons than Cu), and modifies the local structure and chemistry around the Ce atoms (since Sn occupies preferentially the special Cu site which is next to each of the four Ce atoms in the unit cell). The result of these three effects is the progressive dominance of the RKKY interaction over the Kondo interaction and thus the establishment of a long-range antiferromagnetic order of the uncompensated moments. This order is most robust in the ordered material CeCu 5 Sn and therefore higher magnetic fields are required to suppress the magnetic transition.

While further experiments on the crystal/magnetic structure and on the physical properties as well as band structure calculations would be valuable to establish an accurate Doniach's phase diagram of T vs x of CeCu 6-x Sn x , the present results unravelled a very interesting structural chemistry behind the physics of this system and a remarkable interplay between its composition, structure and magnetism. Whether this interplay can be further tuned by control parameters like pressure or higher magnetic fields remains an open interesting question.

x a (Å) b (Å) c (Å) V (Å 3 ) 0 8.

  where I k and I ck are the observed and the calculated integrated nuclear intensity at the k th reflection respectively. Guidelines of the Rietveld refinement can be found elsewhere (ref. 13). The neutron scattering lengths used were b Ce = 0.484 10 -14 m, b Sn = 0.6225 10 -14 m and b Cu = 0.7718 10 - 14 m, values taken from reference 14.

  3a. Crystal structureThe existence of the series RCu 5 Sn (R = La, Ce, Pr, Nd and Sm) has been first reported by and independently byFornasini et al. (ref. 16). According to Skolozdra, the compounds RCu 5 Sn adopt the CeCu 6 type-structure (ref. 17) but without forming a complete solid solution between RCu 6 and RCu 5 Sn. Conventional powder X-ray diffraction of the CeCu 6-x Sn x (x = 0, 0.25, 0.65, 0.75, 0.85 and 1) compounds investigated here, revealed only traces of elemental Sn in the samples with higher Sn content and only traces of CeO 2 in the CeCu 6 sample.The CeCu 6-x Sn x compounds are found to crystallize in the orthorhombic symmetry adopting the Pnma space group similarly to the CeCu 5 Au compound (ref.18).The lattice parameters deduced from the refinement of the neutron diffraction data of the CeCu 6-x Sn x (x = 0, 0.25, 0.65. 0.75, 0.85, 1) compounds are summarized in TableI. As can be seen from this table, the lattice expands with the substitution of Cu by Sn, which indicates that the substitution is equivalent to negative chemical pressure on the lattice. The composition x dependencies of the lattice parameters a, b and c and of the volume V are given in Figure I. The b lattice parameter is almost insensitive to the substitution of Cu by Sn in the CeCu 6-x Sn x structure. On the contrary, both lattice parameters a and c increase continuously with x; c being by far the most sensitive to the Sn concentration. We now focus on the crystal structure of the compounds CeCu 5.25 Sn 0.75 and CeCu 5 Sn which are representative phases in the area of high Sn content, where antiferromagnetic order at low temperatures has been reported (ref. 10). Their crystal structure was investigated in detail above their corresponding Néel temperatures in order to determine the exact site occupancy of Sn in the structure and its influence to the other structural parameters and to compare it with the structure of the parent compound CeCu 6 . An example of the Rietveld refinement fit is given in Figure II for CeCu 5.25 Sn 0.75 at 5 K. The structural results of the Rietveld refinement of the neutron diffraction data are summarized in Tables II and III for the CeCu 5.25 Sn 0.75 and CeCu 5 Sn respectively. It is unambiguously deduced from the refinements that the Sn atoms exhibit a strong preference for the Cu 4c site of the CeCu 6 structure. Such preference for the 4c site confirms earlier reported results

  . The metallic atomic radii used are the following: r Cu = 1.28 Å, r Sn = 1.60 Å, r Ce = 1.846 Å according to the calculations of reference 23. The available atomic cell volume calculated for each atomic site is summarized in Table III for CeCu 5 Sn. It is clear that the position preferentially occupied by Sn atoms is the one exhibiting the larger volume, thus favouring the accommodation of large atoms such as Sn or Au. Similar calculations on the pure CeCu 6 compounds confirm that the position preferred by Sn exhibits the larger W.S. atomic cell volume: 13.7Å 3 against only 12.4 Å 3 and about 12 Å 3 for the Cu_8d and the two other Cu_4c site respectively. Upon substituting Sn for Cu, both the expansion of the unit cell as well as the evolution of the coordinates of the atomic positions permit to accommodate the large Sn atom in the CeCu 6 type crystal structure. Indeed, all atoms of CeCu 6-x Sn x are in general positions of the Pnma space group and therefore the coordinates of the atomic positions are allowed to change. Thus, for the end-member CeCu 5 Sn the volume available for the Sn site reaches 20.79 Å 3 .

  Fig. I: Composition x dependence of the lattice parameters a, b and c and of the unit cell volume V

Fig. II: 5 Sn

 5 Fig. II: Neutron powder diffraction pattern of CeCu 5.25 Sn 0.75 at 5 K. The dots and the full curve refer

2 )

 2 Fig. IV. Specific heat of polycrystalline CeCu 6 at zero magnetic field.

Fig. V :

 V Fig. V: Low temperature specific heat of the CeCu 6-x Sn x compounds with x = 0.25 and x = 0 at zero

Fig

  Fig. VI: Comparison of the specific heat of polycrystalline CeCu 6-x Sn x with high Sn content at zero

  Fig. VII: Evolution of C p /T versus T of the CeCu 5.35 Sn 0.65 for applied magnetic fields of 0, 1, 5, 7

19

 19

  Fig. VIII: Evolution of C p /T versus T of the CeCu 5.25 Sn 0.75 for applied magnetic fields of 0, 3, 5, 7

  Fig. IX: Evolution of C p /T versus T of the CeCu 5 Sn for applied magnetic fields of 0, 3, 5, 7 and 10

20

 20

Table I :

 I Lattice parameters of the CeCu 6-x Sn x compounds at RT.

			117(1)	5.106 (7)	10.182 (2)	422.0
	0.25		8.218 (1)	5.079 (2)	10.281 (1)	429.2
	0.65		8.330 (1)	5.067 (1)	10.551 (2)	445.4
	0.75		8.319 (2)	5.073 (2)	10.520 (3)	444
	0.85		8.328 (1)	5.069 (1)	10.595 (1)	447.2
	1		8.323 (1)	5.064 (1)	10.597 (1)	447.6
	1.2 (ref. 15)	8.342 (1)	5.067 (5)	10.651 (10)	450.2
	Atom	site x		y	z	Occupancy
	Ce	4c 0.230 (3)	0.25	0.561 (2)	1
	Cu1	8d 0.0685 (9) 0.4976 (22) 0.3078 (8) 1
	Cu2	4c 0.0513 (18) 0.25	0.1059 (16) 1
	Cu3	4c 0.3105 (20) 0.25	0.2466 (12) 1
	Cu4	4c 0.4137 (20) 0.25	0.0213 (13) 1
	Sn5/Cu5 4c 0.1416 (19) 0.25	0.8660 (15) 0.70/0.30(7)

Table II :

 II Structural parameters of CeCu 5.25 Sn 0.75 deduced from the refinement of the neutron

powder diffraction data collected at 5 K. The standard deviations are given in parenthesis. Space group: Pnma; lattice parameters: a = 8.297(1) Å, b = 5.067(1) Å, c = 10.555(2) Å. Rietveld refinement agreement factors: R p = 7.81%, R wp = 7.60%, R B = 2.96%.

Acknowledgements

The work of J. P. and P. J. is a part of the research plan MSM 0021620834 that is financed by the Ministry of Education of the Czech Republic. E.G.M. and O.I. thank GSRT and Égide respectively for financial support through the French-Greek cooperation program PLATON #11202UA and ILL for beam time at the instrument D1B (CRG).

The nearest neighbours and the Wigner-Seitz atomic cell volume are given for each atomic position.


F R Petruzielo

Julien Toulouse
email: julien.toulouse@upmc.fr

C J Umrigar
email: cyrusumrigar@cornell.edu

Compact and Flexible Basis Functions for Quantum Monte Carlo Calculations

Molecular calculations in quantum Monte Carlo frequently employ a mixed basis consisting of contracted and primitive Gaussian functions. While standard basis sets of varying size and accuracy are available in the literature, we demonstrate that reoptimizing the primitive function exponents within quantum Monte Carlo yields more compact basis sets for a given accuracy. Particularly large gains are achieved for highly excited states. For calculations requiring non-diverging pseudopotentials, we introduce Gauss-Slater basis functions that behave as Gaussians at short distances and Slaters at long distances. These basis functions further improve the energy and fluctuations of the local energy for a given basis size. Gains achieved by exponent optimization and Gauss-Slater basis use are exemplified by calculations for the ground state of carbon, the lowest lying excited states of carbon with 5 S o , 3 P o , 1 D o , 3 F o symmetries, carbon dimer, and naphthalene. Basis size reduction enables quantum Monte Carlo treatment of larger molecules at high accuracy.

I. INTRODUCTION

In traditional quantum chemistry (QC) calculations, molecular orbitals are often expanded in a combination of contracted Gaussian basis functions and primitive Gaussian basis functions. For each occupied orbital, a contracted function is constructed to reproduce the corresponding atomic orbital from an effectively singleelectron theory such as Hartree-Fock (HF) [1,2], or the natural orbital from a post-HF method [3,4].

While a single primitive Gaussian has incorrect longrange asymptotic behavior, a contracted basis function can reproduce the correct asymptotics over a reasonable range. However, even contracted functions are unable to produce the correct electron-nucleus cusps [5] since they have zero gradient at the origin. Despite these shortcomings, Gaussians are used in QC calculations because they permit analytical evaluation of the two-electron integrals [6].

In contrast to traditional QC methods, quantum Monte Carlo (QMC) calculations [7] enjoy greater wavefunction flexibility by using Monte Carlo integration to evaluate matrix elements. In particular, bases need not be restricted to Gaussians. For calculations employing a potential that diverges at the nucleus, Slater basis functions can exactly reproduce the correct electron-nucleus cusp and long-range asymptotic behavior of the orbitals. In fact, for all-electron QMC calculations, highly accu-rate results have been obtained by employing compact basis sets consisting of Slater functions with optimized exponents [8,9].

Conversely, the basis sets used for non-divergent pseudopotential calculations in QMC have deviated little from typical QC basis sets. For these pseudopotentials, orbitals have no electron-nucleus cusp. In this case, Gaussian basis functions are more appropriate than Slater functions at small electron-nuclear distances but still have incorrect long-range asymptotics.

Contracted and primitive Gaussian functions are frequently splined on a radial grid for QMC. Splining contracted Gaussians presents a definite computational advantage since evaluating polynomials is much cheaper than evaluating a linear combination of Gaussians. In contrast, splining primitive Gaussians provides minimal benefit at best.

We propose two ideas for improving basis sets for pseudopotential calculations in QMC. First, primitive basis function exponents are optimized for each system. This provides greater accuracy with a compact basis for a wide range of chemical environments and excitation levels. To facilitate optimization, the primitive basis functions remain analytic while the contracted functions are splined.

Second, we propose a novel form of primitive basis function appropriate for calculations involving nondiverging pseudopotentials. These primitives, which we call Gauss-Slater (GS) functions, have the short-range behavior of a Gaussian function and the long-range behavior of a Slater function.

The utility of our improvements is demonstrated by calculations for carbon, the lowest lying excited states of carbon with 5 S o , 3 P o , 1 D o , 3 F o symmetries, carbon dimer, and naphthalene. This paper is organized as follows. In Section II, the form and properties of Gauss-Slater functions are introduced. In Section III, results of our calculations are discussed. In Section IV, concluding remarks are provided. In the Appendices, technical details are discussed.

II. GAUSS-SLATER BASIS FUNCTIONS

We define Gauss-Slater (GS) functions as

ϕ ζ nlm (r, θ, φ) = N ζ n r n-1 e -(ζr) 2 1+ζr Z m l (θ, φ), (1)
where r, θ, φ are the standard spherical coordinates, n is the principal quantum number, l is the azimuthal quantum number, m is the magnetic quantum number, N ζ n is the normalization factor, and Z m l (θ, φ) is a real spherical harmonic.

Notice that for r ≪ 1 the GS behaves like a Gaussian:

ϕ ζ nlm (r, θ, φ) ∼ = N ζ n r n-1 e -(ζr) 2 Z m l (θ, φ), (2)
and for r ≫ 1 the GS behaves like a Slater:

ϕ ζ nlm (r, θ, φ) ∼ = N ζ n r n-1 e -ζr Z m l (θ, φ). (3)
The GS drift velocity and local energy are well behaved at long distances, while for Gaussians these quantities diverge. Unlike Gaussians and Slaters, normalization of GSs has no closed form expression. Nevertheless, normalizing an arbitrary GS is trivial with the following scaling relation (see Appendix A) between N ζ n and N 1 n ,

N ζ n = ζ n+1/2 N 1 n . (4)
Since GSs are not analytically integrable, the exponential part must be expanded in Gaussians for use in quantum chemistry programs that employ analytic integrals for evaluating the matrix elements. This expansion is

N ζ n e -(ζr) 2 1+ζr = i c ζ i 2(2α ζ i ) n+ 1 2 Γ(n + 1 2 ) e -α ζ i r 2 , (5)
where c ζ i is the i th expansion coefficient, and α ζ i is the i th Gaussian exponent. Additionally, the following scaling relations (see Appendix A) hold for the expansion coefficients and Gaussian exponents:

α ζ i = ζ 2 α 1 i ( 6
)
c ζ i = c 1 i . (7)
Once the Gaussian expansions are found for unit exponents, expansions of arbitrary GSs follow immediately from the scaling relations.

III. RESULTS

For all applications discussed in this paper, variational Monte Carlo (VMC) and diffusion Monte Carlo (DMC) [15] calculations are performed with the CHAMP QMC code [START_REF] Champ | a quantum Monte Carlo program written[END_REF] and employ the pseudopotentials and accompanying basis sets of Burkatzki, Filippi and Dolg (BFD) [START_REF] Burkatzki | [END_REF]. We choose these pseudopotentials and basis functions since they were constructed for use in QMC and have proved to be quite accurate.

The wavefunction is of the standard Slater-Jastrow form. All wavefunction parameters including Jastrow parameters, Configuration State Function (CSF) coefficients (where applicable), orbital coefficients, and primitive exponents (where applicable) are optimized via the linear method [12][13][14]. Optimization is performed on a linear combination of the energy and variance of the local energy with weights 0.95 and 0.05, respectively. Optimizing just the energy yields slightly lower energies and somewhat higher variances.

For each system considered, calculations are performed with three different basis sets: (1) the BFD basis, (2) fixed contracted functions and analytical Gaussian primitives with optimized exponents, and (3) fixed contracted functions and analytical Gauss-Slater primitives with optimized exponents. We refer to these basis sets as BFD, G, and GS, respectively.

These three cases allow us to evaluate the improvements our two methods provide to the current basis sets used in QMC. First, if both the G and GS basis sets significantly outperform the BFD basis, then the utility of reoptimizing the basis exponents within QMC will be established. Second, the utility of the GS basis depends on its performance relative to the G basis.

A. Ground State Carbon Atom

For the carbon atom ground state, 3 P , we consider a complete active space (CAS) wavefunction with an active space generated by distributing the four valence electrons among the thirteen orbitals of the n = 2 and n = 3 shells. Denoted by CAS (4,13), this wavefunction consists of 83 CSFs comprised of 422 determinants.

In general, a single Slater determinant will not be a CSF when a certain number of electrons have been excited relative to the ground state HF Slater determinant. However, a CSF can be produced from an arbitrary Slater determinant by applying projection operators for angular momentum L and spin Ŝ. Since states with the same L and S but different L z and S z are degenerate, we are free to choose convenient L z and S z states. We choose L z = 0 to make the wavefunctions real to within a position independent phase, and we choose S z = S to yield the minimum number of determinants in the CSF. Since the carbon ground state has L = 1, S = 1, the projection operators are of the form

PL = L ′ =1 L -L ′ (L ′ + 1) (8) PS = S ′ =1 Ŝ -S ′ (S ′ + 1) , (9)
where the product over all possible angular momentum and spin values omits the desired L = 1 and S = 1 values. Carbon atom VMC results for each basis set are shown in Table I. Included for comparison, coupled cluster calculations with single and double excitations and perturbative triple excitations (CCSD(T)) values for the BFD basis [START_REF] Burkatzki | [END_REF] exhibit much larger dependence on basis size than QMC results.

Both the G and GS basis sets outperform the BFD basis set. The 2z G basis exhibits a modest gain of 0.3 mH in energy compared to the corresponding BFD basis. The 2z GS basis exhibits larger gains of 1 mH in energy and 28 mH in σ, the root-mean-square (RMS) fluctuations of the local energy. The 3z GS basis yields identical results, within statistical error, to the 5z BFD basis.

Carbon atom DMC results for each basis set are shown in Table II. These calculations were performed with a time step of τ = 0.01 H -1 which leads to negligible time step error for these high quality wavefunctions. DMC depends less on basis size than VMC, as is immediately apparent from the data. Nevertheless, both the G and GS basis sets outperform the BFD basis set. The 3z GS basis yields identical results, within statistical error, to the 5z BFD basis, and an energy 0.1 mH lower than the 3z BFD basis.

Both VMC and DMC results indicate that reoptimizing primitive basis function exponents provides improvements which can be significant for the GS basis. In large systems, the ability to use a 3z basis in place of a 4z or 5z basis determines whether a calculation can be performed.

B. Excited States of Carbon

We consider the lowest lying excited states of carbon with 5 S o , 3 P o , 1 D o , and 3 F o symmetries. These states have configurations 2s 1 2p 3 , 2s 2 2p 1 3s 1 , 2s 2 2p 1 3d 1 , and 2s 2 2p 1 3d 1 , respectively. The 3 P o , 1 D o , and 3 F o states have much higher energy than the ground state and 5 S o excited state.

The dominant CSF for each of these three states occupies orbitals that are unoccupied in the HF ground state. For fair comparison, the BFD basis therefore must be augmented. The diffuse functions of the aug-cc-pVnZ basis sets [4,16,17] are used for this purpose. The BFD nz basis then becomes an (n + 1)z basis.

Obtained by application of the projection operators discussed in Section III A, the dominant CSF for each of the four excited states has one, one, four, and three Slater determinants, respectively. VMC results for energies and σ of each system are shown in Figures 1 and2. In all cases, the reoptimized exponents provide significant gains in both energy and σ. Results for the three higher lying states demonstrate that reoptimized exponents are essential for describing states containing orbitals unoccupied in the HF ground state. For these systems, 2z results using the G and GS basis sets are substantially better than 5z results for the augmented BFD basis set. In the most extreme case of 3 F o , the 2z GS basis results in 30 mH lower energy and 110 mH lower σ than the 5z BFD basis.

The importance of the reoptimized exponents is evident for the excited states of carbon. However, benefits of the GS basis relative to the G basis are never more than several tenths of a mH. On the scale of the plots in Figures 1 and2, many G and GS basis results coincide. For each n, the nz basis consists of n S functions, n P functions, and n -1 D functions (where applicable). For 3 P o , 1 D o , 3 F o calculations, the BFD basis is augmented with diffuse functions of the aug-cc-pVnZ basis sets [4,16,17]. For G and GS basis sets, only 2z and 3z calculations were performed. In many cases, results for G and GS bases are indistinguishable on this scale.

C. Carbon Dimer

Single determinant calculations were performed for the carbon dimer with initial wavefunctions generated from the QC code GAMESS [18].

VMC results for each basis set are shown in Table III. The G and GS basis sets outperform the BFD basis set. In particular, the 2z G basis attains a 0.6 mH lower energy than the corresponding BFD basis, and the GS basis yields a 3.2 mH lower energy than the BFD basis. The 3z GS basis yields an energy within 0.3 mH of and a σ identical to the 5z BFD basis results.

DMC results for each basis set are shown in Table IV.

These calculations were performed with a time step of τ = 0.005 H -1 which leads to negligible time step error. The 2z G and GS basis sets significantly outperform the corresponding BFD basis set. The 2z GS basis yields a result that is essentially converged with respect to basis size.

Both VMC and DMC results indicate that reoptimizing primitive basis function exponents provides improvements which can be significant for the GS basis.

D. Naphthalene

Single determinant calculations were performed for naphthalene, C 10 H 8 , with initial wavefunctions generated from the QC code GAMESS [18].

Calculations were performed only for the 2z basis. The intention of this section is not to produce an energy converged with respect to basis size, but to demonstrate that the utility of reoptimizing primitive basis functions, and GS primitives in particular, extends to large systems. VMC and DMC results for each basis set are shown in Table V. The DMC calculations were performed with a time step of τ = 0.01 H -1 .

At both the VMC and DMC level, wavefunction quality increases by reoptimizing the primitive Gaussian exponents. The GS basis provides further improvement. In particular, even at the DMC level, the 2z G basis attains a 4 mH lower energy than the corresponding BFD basis, and the GS basis yields a 15 mH lower energy than the BFD basis. This is significant since DMC is less sensitive to basis set choice than VMC. First, we reoptimized the primitive basis function exponents for each system because the exponents of standard QC and QMC basis sets, such as the Burkatzki, Filippi and Dolg (BFD) basis, represent a compromise.

These standard exponents are designed to yield good energies for some range of chemical environments and excitation levels, but they cannot be optimal for all systems. We have shown that reoptimizing primitive basis function exponents for each system yields significant improvements in the energy and fluctuations of the local energy, σ. The most pronounced benefits were observed in higher-lying excited state calculations. In the most extreme case of 3 F o at the variational Monte Carlo (VMC) level, the 2z mixed basis was 30 mH lower in energy and 110 mH lower in σ than the 5z numerical basis. Although not discussed in this paper, we have found that reoptimization of standard Slater basis exponents used in all-electron calculations also provides considerable improvements in energy and σ.

Second, we introduced Gauss-Slater (GS) basis functions for non-divergent pseudopotential calculations. GS functions behave like Gaussians at short distances and Slaters at long distances. In all systems considered, results obtained using a mixed basis comprised of contracted and primitive basis improved when optimized Gaussian primitives were replaced by optimized GS primitives. Importantly, for carbon dimer at the DMC level the 2z GS total energies are nearly converged with respect to basis size.

A 3z mixed basis with optimized GSs for carbon atom or carbon dimer produces results comparable to the 5z BFD basis. Since the number of orbital coefficients to be optimized scales quadratically with basis size, the use of a more compact basis allows larger problems to be attacked in QMC.

FIG. 1 :

 1 FIG.1: VMC energies in Hartrees for lowest lying excited states of carbon with 5 S o , 3 P o , 1 D o , 3 F o symmetries. For each n, the nz basis consists of n S functions, n P functions, and n -1 D functions (where applicable). For 3 P o , 1 D o , 3 F o calculations, the BFD basis is augmented with diffuse functions of the aug-cc-pVnZ basis sets[4,16,17]. For G and GS basis sets, only 2z and 3z calculations were performed. In many cases, results for G and GS bases are indistinguishable on this scale.

FIG. 2 :

 2 FIG. 2: RMS fluctuations of VMC local energies. See Figure 1 for notation and details.

TABLE I :

 I VMC energy and RMS fluctuations of the local energy, σ, in Hartrees for CAS(4,13) ground state of carbon using BFD, G, and GS basis functions. Statistical errors on the last digit are shown in parentheses. For each n, the nz basis consists of n S functions, n P functions, and n -1 D functions. CCSD(T) values for the BFD basis are included for comparison[START_REF] Burkatzki | [END_REF].

	Type	Size Energy (H) σ (H)
	BFD	2z -5.43161(3) 0.1395(6)
		3z -5.43306(2) 0.099(3)
		4z -5.43332(2) 0.0904(2)
		5z -5.43341(2) 0.0905(4)
	G	2z -5.43196(3) 0.138(2)
		3z -5.43324(2) 0.0989(5)
	GS	2z -5.43264(2) 0.1114(4)
		3z -5.43344(2) 0.0898(2)
	CCSD(T) 2z -5.409230 N/A
		3z -5.427351 N/A
		4z -5.431486 N/A
		5z -5.432494 N/A
	Type Size Energy (H)
	BFD 2z -5.43314(2)
		3z -5.43395(2)
		4z -5.43404(1)
		5z -5.43407(1)
	G	2z -5.43342(2)
		3z -5.43400(2)
	GS 2z -5.43356(2)
		3z -5.43407(1)

TABLE II: DMC energy in Hartrees for CAS

(4,13)

ground state of carbon using BFD, G, and GS basis functions. Statistical errors on the last digit are shown in parentheses. For each n, the nz basis consists of n S functions, n P functions, and n -1 D functions. Calculations were performed with a time step of τ = 0.01 H -1 which leads to a negligible time step error for these high quality wavefunctions.

TABLE III :

 III VMC energy and σ in Hartrees for single determinant carbon dimer ground state. Statistical errors on the last digit are shown in parentheses. For each n, the nz basis consists of n S functions, n P functions, n -1 D functions, n -2 F functions, and n -3 G functions.

	Type Size Energy (H) σ (H)
	BFD 2z -11.02644(4) 0.4343(9)
		3z -11.03003(4) 0.4172(4)
		4z -11.03094(4) 0.4127(7)
		5z -11.03095(4) 0.4113(6)
	G	2z -11.02707(4) 0.4288(7)
		3z -11.03030(4) 0.4183(6)
	GS 2z -11.02968(4) 0.4191(6)
		3z -11.03065(4) 0.4109(6)

TABLE IV

 IV

		005 H -1
	which leads to negligible time step error.
	Type Size Energy (H)
	BFD 2z -11.05561(3)
		3z -11.05719(4)
		4z -11.05728(4)
		5z -11.05723(4)
	G	2z -11.05632(3)
		3z -11.05717(3)
	GS 2z -11.05702(4)
		3z -11.05719(3)

: DMC energy in Hartrees for ground state of carbon dimer using BFD, G, and GS basis functions. Statistical errors on the last digit are shown in parentheses. For each n, the nz basis consists of n S functions, n P functions, n -1 D functions, n -2 F functions, and n -3 G functions. Calculations were performed with a time step of τ = 0.

TABLE V :

 V VMC energy and σ, and DMC energy in Hartrees for single determinant ground state naphthalene, C10H8. Statistical errors on the last digit are shown in parentheses. DMC calculations were performed with a time step of τ = 0.01 H -1 . For carbon, the 2z basis includes 2 S function, 2 P function, 1 D function. For hydrogen, the 2z basis includes 2 S functions and 1 P function. Calculations for the ground state of carbon, the lowest lying 5 S o , 3 P o , 1 D o , 3 F o excited states of carbon, carbon dimer, and naphthalene demonstrate the utility of our contribution.

	Type Size VMC Energy (H) VMC σ (H) DMC Energy (H)
	BFD 2z -61.5193(5)	0.980(1)	-61.6479(5)
	G	2z -61.5273(4)	0.938(1)	-61.6518(5)
	GS 2z -61.5438(4)	0.927(2)	-61.6634(5)
		IV. CONCLUSION
		Quantum Monte Carlo (QMC) methods have the ad-
	vantage over standard quantum chemistry (QC) methods
	of rapid convergence with increasing basis size. Basis-
	size dependence for pseudopotential calculations in QMC
	is further reduced by two basis set improvements intro-
	duced in this work.

TABLE VIII :

 VIII Basis exponents for the lowest lying excited states of carbon with 5 S o , 3 P o , 1 D o , 3 F o symmetries. For each n, the nz basis includes n -1 S primitives, n -1 P primitives, and n -1 D primitives (where applicable).Type Size L 5 S o Exp. 3 P o Exp. 1 D o Exp. 3 F o Exp.

	G	2z S 0.110	0.006	0.087	0.112
		P 0.150	0.157	0.109	0.766
		D N/A	N/A	0.006	0.006
		3z S 0.356	0.010	0.094	0.096
		S 2.145	0.284	0.496	0.689
		P 0.809	4.868	0.079	0.117
		P 2.262	6.707	0.854	0.524
		D N/A	N/A	0.007	0.007
		D N/A	N/A	0.407	0.822
	GS 2z S 1.960	0.059	0.756	0.570
		P 1.860	3.002	0.763	0.588
		D N/A	N/A	0.351	0.293
		3z S 1.148	0.371	0.160	0.894
		S 1.451	0.376	0.931	2.045
		P 1.137	0.691	0.043	1.326
		P 1.325	2.087	0.444	3.016
		D N/A	N/A	0.238	0.295
		D N/A	N/A	0.953	1.545

TABLE IX :

 IX Basis exponents for ground state of naphthalene, C10H8, using G and GS basis functions. For carbon, the 2z basis includes 1 S primitive, 1 P primitive, 1 D primitive. For hydrogen, the 2z basis includes 1 S primitive.

	Type Size L C10H8 C Exp. C10H8 H Exp.
	G	2z S 0.139	0.099
		P 0.191	N/A
		D 0.754	N/A
	GS 2z S 0.875	0.798
		P 1.118	N/A
		D 2.109	N/A

V. ACKNOWLEDGMENTS

This work was supported by the NSF (grant DMR-0908653) and by the DOE (grant DOE-DE-FG05-08OR23336). Computations were performed in part at the Cornell NanoScale Facility, a member of the National Nanotechnology Infrastructure Network, and, at the Computation Center for Nanotechnology Innovation at Rensselaer Polytechnic Institute.

(A3)

Hence, the scaling relation for the normalization factor is

The values of N 1 n are given in Table VI. To derive the scaling relations for the parameters α ζ i and c ζ i in the Gaussian expansion of the Gauss-Slater functions, suppose the best-fit expansion for ζ = 1 is

Using Eqn. (A4) and performing the substitution r → ζr results in

where

Appendix B: Spatial Derivatives A general unnormalized radial basis function has the form

where g ζ (r) is an arbitrary function. The gradient is

where

The Laplacian is

where

For Gauss-Slater functions,

For Gaussian functions,

For Slater functions,

Appendix C: Parameter Derivatives Wavefunction optimization via the linear method requires both the derivatives of the wavefunction with respect to the exponent parameters ζ, and the Hamiltonian acting on those derivatives. From Eqn. (A4), the derivative of the normalization with respect to the exponent is

Now consider a general unnormalized radial basis function of the form

where g ζ (r) is an arbitrary function. The derivative of the radial part of the wavefunction with respect to the exponent is

where

The gradient is

The Laplacian is

r). (C6)

For Gauss-Slater functions,

For Gaussian functions,

For Slater functions,

Appendix D: Exponents

To promote use of this basis, exponents for each system studied are given in this appendix. Only exponents of the primitives are given, as the contracted functions are presented elsewhere [START_REF] Burkatzki | [END_REF]. Exponents for the ground state of carbon using CAS(4,13) CSFs, and, the ground state of carbon dimer using a single CSF are shown in Table VII. Exponents for the lowest lying excited states of IX. In naphthalene, since atoms of the same atomic species are located at inequivalent geometrical locations, one could independently optimize the exponents for each inequivalent atom, but we have not done so because we expect the resulting gain to be small. We have found that the carbon S and P exponents change relatively little from one molecule to another (though they do differ more for the atom) while there is considerable leeway in the D exponents (they change considerably even from one optimization to another for a given molecule). This is because the energy and σ are not as sensitive to the D basis functions as they do not appear in the ground-state determinant of the carbon atom. Hence it is possible to find an approximately optimal set of exponents for the atoms in a large molecule, by optimizing them for a small molecule with the same atoms.