François Bottin
email: francois.bottin@cea.fr

Jordan Bieder

Johann Bouchet

a-TDEP Temperature Dependent Effective Potential for Abinit -Lattice dynamic properties including anharmonicity

Keywords: Phonons, Anharmonicity, ab initio, Molecular Dynamics, Thermodynamics

In this paper, we present the a-TDEP post-process code implemented in the Abinit package. This one is able to capture the explicit thermal effects in solid state physics and to produce a large number of temperature dependent thermodynamic quantities, including the so-called anharmonic effects. Its use is straightforward and require only a single ab initio molecular dynamic (AIMD) trajectory. A Graphical User Interface (GUI) is also available, making the use even easier.

. In particular, we present the various algorithms and schemes used in a-TDEP which enable to obtain the effective Interatomic Force Constants (IFC). The 2 nd and 3 rd order effective IFC are produced self-consistently using a leastsquare method, fitting the AIMD forces on a model Hamiltonian function of the displacements. In addition, we stress that we face to a constrained leastsquare problem since all the effective IFC have to fulfill the several symmetry rules imposed by the space group, by the translation or rotation invariances of the system and by others.

Numerous thermodynamic quantities can be computed starting from the

Introduction

The study of thermal properties of materials has always been the subject of a significant research effort in solid state physics. Whatever the purpose, experiments and theory have brought together their discoveries to understand the mechanisms that governed temperature effects. A large number of fields of research are concerned and interested in a still better description of materials when they undergo some temperature variations : astrophysics, superconductivity, aerospace industry, electronics, geoscience... All of them require the building of phase diagrams, the refinement of Equation of States (EoS) and the characterization of thermodynamic quantities.

The theory of lattice dynamics has been developed for this purpose and specifies that thermal effects can be captured by describing precisely the lattice vibrations of the system. In this framework, the key quantity is the phonon frequencies ω(V, T ) which depend on the volume V and the temperature T . Within the harmonic approximation (HA), the most famous approach based on the quantum harmonic crystal, the lattice vibrations ω HA (V ) are temperature independent and the temperature effects are only taken into account through the population of phonon states by the Bose-Einstein statistics. If this approximation has encountered many success, nevertheless this one failed to describe the thermal expansion of materials. That is the reason why another approximation has been proposed : the quasi-harmonic approximation (QHA). This one considers that the phonon frequencies depend on the temperature through the variation of the volume : ω QHA (V (T )). Even if the temperature is not explicitly taken into account, the QHA is able to capture most of standard thermal effects encountered by materials.

However, the QHA reaches its limits when the atoms explore the anharmonic part of the potential. If this can happen at low temperature (for hydrogen and helium, for instance), this is generally encountered when materials are strongly heated, when some crystalline phases are (de)stabilized wrt temperature, when thermodynamic conditions are close to phase transitions... In these cases, an harmonic description of the potential is no longer sufficient. Moreover, if the QHA already includes some effects going beyond the harmonic approximations (the thermal expansion, for instance), it is not able to account for thermal effects at constant volume. In this case, the temperature effects not described by the harmonic approximation can be split in two parts [START_REF] Anderson | The grüneisen ratio for the last 30 years[END_REF][START_REF] Lucazeau | Effect of pressure and temperature on raman spectra of solids: anharmonicity[END_REF][START_REF] Holzapfel | Effects of intrinsic anharmonicity in the mie-grüneisen equation of state and higher order corrections[END_REF][START_REF] Ambrumenil | Lattice dynamics and negative thermal expansion in the framework compound ZnNi(CN) 4 with two-dimensional and three-dimensional local environments[END_REF] :

• one described by the QHA, named isothermal or extrinsic, able to describe the thermal expansion but inefficient to depict the thermal behaviour of materials at constant volume,

• another going beyond the QHA, named isochoric or intrinsic, which treats explicitly the dependency wrt the temperature and reproduces thermal effects at constant volume : ω Anh (V (T ), T ).

One can formalize and synthesize this statement as follows :

∂ω ∂T p = ∂ω ∂T V + ∂ω ∂V T ∂V ∂T p (1)
The second term of the right hand side of the previous equation is included in the QHA, but not the first one only treated when the temperature dependency is explicit.

As well for material science as for theoretical issues, the calculation of lattice vibrations (whatever the thermodynamic conditions) is becoming one of the most important challenge in solid state physics. For this purpose, the ab initio calculations are well suited : both the effects of pressure and temperature can be taken into account, without any assumption about the interatomic potential. Two approaches have been developed in order to compute the phonon frequencies at the harmonic level : either by brute force using finite differences [START_REF] Parlinski | First-principle determination of the soft mode in cubic zro 2[END_REF][START_REF] Togo | First-principles calculations of the ferroelastic transition between rutile-type and cacl 2 -type sio 2 at high pressure[END_REF] or by employing a more elegant method, the density functional perturbation theory (DFPT) [START_REF] Gonze | Dynamical matrices, born effectives charges, dielectric permittivity tensors, and interatomic force constants from densityfunctional perturbation theory[END_REF][START_REF] Baroni | Phonons and related crystal properties from density-functional perturbation theory[END_REF]. If recent developments show that finite differences ab initio calculations are able to capture strong anharmonic effects [START_REF] Feng | Quantum mechanical prediction of four-phonon scattering rates and reduced thermal conductivity of solids[END_REF][START_REF] Tadano | First-principles lattice dynamics method for strongly anharmonic crystals[END_REF] it seems more challenging for DFPT to compute the terms above the cubic one [START_REF] Paulatto | Anharmonic properties from a generalized third-order ab initio approach: Theory and applications to graphite and graphene[END_REF] using the 2n + 1 theorem [START_REF] Gonze | Density-functional approach to nonlinearresponse coefficients of solids[END_REF].

Since the beginning of the 60's, several microscopic formalisms has been proposed in order to go beyond the simple quantum harmonic crystal and to deal with the anharmonic effects [START_REF] Leibfried | Theory of anharmonic Effects in Crystals[END_REF][START_REF] Maradudin | Scattering of neutrons by anharmonic crystal[END_REF][START_REF] Cowley | Anharmonics crystals[END_REF][START_REF] Barron | Dynamical properties of solids : Volume 1. Crystalline Solids, fundamentals[END_REF][START_REF] Brüesch | Phonons : Theory and Experiments I[END_REF][START_REF] Wallace | Thermodynamics of Crystals[END_REF][START_REF] Wallace | Statistical physics of crystals and liquids[END_REF]. Not only do they offer a description of lattice vibrations going beyond the QHA (that is to say capturing thermal intrinsic effects) but also they give a comprehensive framework including anharmonicity in thermodynamics and elasticity. If these theoretical developments were very successful in the past sixty years, they were restricted to analytical results. Nowadays, due to the increasing power of supercomputers, they can be introduced in computational codes in order to access quantities unreachable in the past.

A large number of theoretical approaches have been developed during the last ten years in order to capture the anharmonics effects starting from ab initio calculations. These ones have led to several computational codes which have recently emerged in the solid state physics community. Without claiming to make an exhaustive list, we can mention :

• the so-called "Self-Consistent Ab-Initio Lattice Dynamics" (SCAILD) method developed by Souvatzis et al. [START_REF] Souvatzis | Entropy driven stabilization of energetically unstable crystal structures explained from first principles theory[END_REF][START_REF] Souvatzis | Dynamical stabilization of cubic zro 2 by phonon-phonon interactions: Ab initio calculations[END_REF][START_REF] Souvatzis | The selfconsistent ab initio lattice dynamical method[END_REF][START_REF] Luo | Dynamical stability of body center cubic iron at the earths core conditions[END_REF],

• the famous python code for phonons calculations (PHONOPY) developed by Togo and Tanaka [START_REF] Togo | First-principles calculations of the ferroelastic transition between rutile-type and cacl 2 -type sio 2 at high pressure[END_REF][START_REF] Togo | First principles phonon calculations in materials science[END_REF][START_REF] Linnera | Ab initio study of the lattice thermal conductivity of cu 2 o using the generalized gradient approximation and hybrid functional methods[END_REF] and the new code (DynaPhoPy) able to extract phonon quasiparticles from AIMD simulations [START_REF] Carreras | Dynaphopy : A code for extracting phonon quasiparticles from molecular dynamics simulations[END_REF],

• the computational code named "Anharmonic LAttice MODEl" (ALAM-ODE) developed by Tadano et al. [START_REF] Tadano | Anharmonic force constants extracted from first-principles molecular dynamics: applications to heat transfer simulations[END_REF][START_REF] Tadano | Self-consistent phonon calculations of lattice dynamical properties in cubic srtio 3 with first-principles anharmonic force constants[END_REF][START_REF] Tadano | Impact of rattlers on thermal conductivity of a thermoelectric clathrate: A first-principles study[END_REF][START_REF] Tadano | First-principles lattice dynamics method for strongly anharmonic crystals[END_REF] • the "Stochastic Self-Consistent Harmonic Approximation" (SSCHA) developed by Errea et al. [START_REF] Errea | First-principles theory of anharmonicity and the inverse isotope effect in superconducting palladium-hydride compounds[END_REF][START_REF] Errea | Anharmonic free energies and phonon dispersions from the stochastic self-consistent harmonic approximation: Application to platinum and palladium hydrides[END_REF][START_REF] Paulatto | First-principles calculations of phonon frequencies, lifetimes, and spectral functions from weak to strong anharmonicity: The example of palladium hydrides[END_REF][START_REF] Borinaga | Anharmonic effects in atomic hydrogen: Superconductivity and lattice dynamical stability[END_REF][START_REF] Bianco | Second-order structural phase transitions, free energy curvature, and temperaturedependent anharmonic phonons in the self-consistent harmonic approximation : Theory and stochastic implementation[END_REF],

• the approach named "Compressive Sensing Lattice Dynamics" developed by Nelson, Zhouet al. [START_REF] Nelson | Compressive sensing as a paradigm for building physics models[END_REF][START_REF] Zhou | Lattice anharmonicity and thermal conductivity from compressive sensing of first-principles calculations[END_REF],

• the method proposed by Glensk et al. allowing to deduce the anharmonic contributions via a direct derivation of the Gibbs energy [START_REF] Glensk | Breakdown of the arrhenius law in describing vacancy formation energies: The importance of local anharmonicity revealed by ab initio thermodynamics[END_REF][START_REF] Glensk | Understanding anharmonicity in fcc materials: From its origin to ab initio strategies beyond the quasiharmonic approximation[END_REF][START_REF] Duff | Improved method of calculating ab initio high-temperature thermodynamic properties with application to zrc[END_REF],

• the Automatic Anharmonic Phonon Library (AAPL) developped by Curtarolo et al. [START_REF] Toher | High-throughpout computational screening of thermal conductivity, debye temperature, and grüneisen parameters using a quasiharmonic debye model[END_REF][START_REF] Plata | An efficient and accurate framework for calculating lattice thermal conductivity of solids : Aflow-appl automatic anharmonic phonon library[END_REF] In this work, we present an implementation of the method originally proposed by Esfarjani and Stokes [START_REF] Esfarjani | Method to extract anharmonic force constants from first principles calculations[END_REF][START_REF] Esfarjani | Heat transport in silicon from first-principles calculations[END_REF][START_REF] Shiomi | Thermal conductivity of half-heusler compounds from first-principles calculations[END_REF] then developed by Hellman and collaborators [START_REF] Hellman | Lattice dynamics of anharmonic solids from first principles[END_REF][START_REF] Hellman | Temperature dependent effective potential method for accurate free energy calculations of solids[END_REF][START_REF] Hellman | Temperature-dependent effective third-order interatomic force constants from first principles[END_REF][START_REF] Shulumba | Lattice thermal conductivity of polyethylene molecular crystals from first-principles including nuclear quantum effects[END_REF]. This one, named "Temperature Dependent Effective Potential" (TDEP), allows to extract the effective interatomic force constants (IFC) by means of AIMD simulations or through a stochastic generation of the configurations. As already highlighted by other authors, this method can be applied successfully to go beyond QHA and capture explicit thermal effects [START_REF] Bouchet | Thermal evolution of vibrational properties of α-u[END_REF][START_REF] Dewaele | High pressure-temperature phase diagram and equation of state of titanium[END_REF][START_REF] Bouchet | High-temperature and high-pressure phase transitions in uranium[END_REF][START_REF] Dorado | Phonon spectra of plutonium at high temperatures[END_REF][START_REF] Dewaele | Iron under conditions close to the α -γ -triple point[END_REF][START_REF] Bouchet | Ab initio calculations of the B1-B2 phase transition in MgO[END_REF].

If a large number of equations and definitions are already published in the literature, few details are available concerning the implementation of this method. In the following, we propose to give them in the particular case of our home-made implementation in the Abinit code [START_REF] Gonze | First-principles computation of material properties: the abinit software project[END_REF]. In the first section, we present the TDEP method as it is exposed in the literature. Then, we detail its implementation in the ABINIT code and focus on the algorithms used in practice. In the following, we describe all the capabilities of the a-TDEP code ; i.e. all the thermodynamic quantities that can be computed. Finally, we show the results obtained for three prototypical materials : Si, Zr and U.

The Temperature Dependent Effective Potential

The Interatomic Forces Constants

Let us define a 3-dimensional crystal composed of N a atoms i in their equilibrium positions which undergo small displacements u α i along the α direction. The potential energy of this system can be rewritten using a Taylor expansion around the equilibrium :

U = U 0 + p≥1 1 p ! α 1 ...αp i 1 ...ip (p) Φ α 1 ...αp i 1 ...ip p k=1 u α k i k (2)
with U 0 the minimum of the potential energy. In this equation, and in the following, the Latin letters in subscripts i, j, k, ... (but also i 1 ,...,i p ) and the Greek letters in superscripts α, β, γ, ... (but also α 1 , ..., α p ) will define the atoms and the cartesian directions, respectively. In Equation 2, the first term of the development

(1)
Φ α i defines the net force acting on the atom i and is null when the expansion is performed around equilibrium positions (corresponding to the thermodynamic conditions).

At the second order, the hamiltonian includes all pair interactions between atoms i and j. This quantity

Φ αβ ij defines the harmonic part of the potential and is called "the second order IFC" (IFC tot 2 ). In the same way, the third and fourth order IFC (IFC tot 3 , IFC tot 4 ), also called cubic and quartic force constants, are

Φ αβγ ijk and (3)
Φ αβγδ ijkl . These ones contain all the three and four body interactions between the i, j, k and l atoms. When the third and fourth order terms become significant, the system acquires a truly anharmonic behavior with an asymmetric shape of the potential [START_REF] Feng | Four-phonon scattering significantly reduces intrinsic thermal conductivity of solids[END_REF][START_REF] Yue | Metric for strong intrinsic fourth-order phonon anharmonicity[END_REF][START_REF] Ravichandran | Unified first-principles theory of thermal properties of insulators[END_REF].

In this framework, the total force acting on each atom i and for each direction α is written as :

F α 1 i 1 = - p≥2 1 (p -1) ! α 2 ...αp i 2 ...ip (p) Φ α 1 α 2 ...αp i 1 i 2 ...ip p k=2 u α k i k (3)

The principles of the TDEP method

In the TDEP method, the aim is to extract the IFC in an effective way. Now, let us assume that the previous system evolves using AIMD, and consider that N t configurations are acquired along the trajectory. At each time step t, we can store the new positions of the atoms τ α M D,i (t), the forces F α M D,i (t) and define the displacements wrt the equilibrium positions

u α M D,i (t) = τ α M D,i (t) -τ α M D,i
, with τ α M D,i the equilibrium positions. Knowing the forces and displacements, and using Equation 3, it is possible to extract the IFC by solving the following system of equations :

F α 1 M D,i 1 (t) = - P p≥2 1 (p -1) ! α 2 ...αp i 2 ...ip (p) Θ α 1 α 2 ...αp i 1 i 2 ...ip p k=2 u α k M D,i k (t) (4)
The expansion is now performed up to the P th order (in practice up to the 2 nd , the 3 rd or the 4 th order). The fitting procedure employed to solve this system of equations modifies the terms below the truncation by including (in an effective way) the anharmonic contributions coming from the terms above the truncation. Therefore, the IFC are no longer constant and become temperature dependent. That is the reason why we change the notation : in the following, the

Φ α...δ i...l will be referred to as the "true IFC" and the

(p) Θ α...δ i...l
as the "effective IFC". By rewriting Eq. 4 as a function of all IFC coefficients (we will show in the following section how to achieve it) we otain :

F α M D,i (t) = P p≥2 λ f α i,pλ (u M D (t))θ pλ (5)
with θ pλ the λ th coefficient of the IFC tot p matrix and f α i,pλ (u M D (t)) the function including all the displacements (but also the factorial and the minus sign).

The system of equations is now linear and has the form A.x = b, with A the function including all the displacements, b the forces and x the IFC coefficients to determine. When the MD trajectory is sufficiently long, this system is overdetermined; i.e. there are more equations than unknowns. One can solve this system of equations by searching its least squares solution. Let us define the residual :

R = F M D -f .Θ (6)
such as

F M D ≡ F α M D,i (t), f ≡ f α i,pλ (u M D (t)
) and Θ ≡ θ pλ . One measure of smallness of R is to choose θ pλ such that the sum of squares of residual S is as small as possible :

S = min(R T .R) = ||F M D -f .Θ|| 2 (7)
The solution giving the lowest residual (i.e. the IFC coefficients giving the best fit of the MD forces) is the following least squares solution :

Θ = f † .F M D (8)
with f † the pseudoinverse of the f matrix. This latter is a generalization of the inverse matrix and we will see at the end of the next section how to compute it.

How many coefficients in the effective IFC?

If the system is composed of N a atoms, the effective IFC contain (3N a ) p coefficients at the p th order, so one has to compute 4 p=1 (3N a ) p coefficients at the fourth order. On the other hand, at each time step of the MD simulation one accumulates 3N a equations. Therefore, if one wants to acquire "at least" as many equations as unknowns, the MD trajectory have to contain N t = 4 p=1 (3N a ) p-1 time steps. For instance, if one considers a supercell with 100 atoms, it is needed to have "at least" N t =90 301 time steps in order to compute the 27 090 300 coefficients up to the third order.

These few lines unveil the problematic. The calculation of the whole IFC tot p matrices is possible but at a prohibitive computational cost. Moreover, the goal is to have an overdetermined system of linear equations in order to solve it by means of a least-squares method.

Here, we show an expansion of the second order effective IFC, solution of the system of equations :

             F x 1 (t) F y 1 (t) F z 1 (t) F x 2 (t) F y 2 (t) F z 2 (t) . . . F z Na (t)              = -              θ 21 0 0 0 θ 21 0 0 0 θ 21     θ 22 θ
                        u x 1 (t) u y 1 (t) u z 1 (t) u x 2 (t) u y 2 (t) u z 2 (t) . . . u z Na (t)             
The whole 3N a × 3N a matrix is symmetric. Each 3 × 3 pair interaction (see for example the matrix) can exhibit some symmetric, antisymmetric, non-symmetric or null terms. Actually, the thousands of coefficients of the IFC tot p can be strongly reduced. If a significant effort is made upstream of the resolution process by taking into account all the symmetries, the method shown previously then becomes tractable. The five symmetry rules, responsible for the decrease of several orders of magnitude, are detailed in the following.

How to reduce the number of IFC coefficients?

The IFC tensors are symmetric

The energy conservation law impose that the IFC tot p tensors comply with the usual transposition rule :

(2)

Θ αβ ij = (2) Θ βα ji (3) Θ αβγ ijk = (3)
Θ βγα jki =

(3)

Θ γαβ kij = (3)
Θ αγβ ikj =

Θ βαγ jik =

(3)

Θ γβα kji (4)
Θ αβγδ ijkl = (4)
Θ βγδα jkli =

Θ γδαβ klij = (4)
Θ αγδβ iklj = (4)
Θ βαγδ jikl =

Θ γδβα klji = ... This invariance under the permutation of the indices is easy to implement. The total number of coefficients is divided by 2 for the IFC tot 2 , 6 for the IFC tot 3 and 24 for the IFC tot 4 .

The acoustic sum rule

The momentum conservation principle implies that there is no mass center displacement. To write this condition properly we split up the atomic index and expand the atomic coordinates as follows : τ i (a) = r i + R(a) with r i the atomic position of the atom i in the cell 0 and R(a) the lattice vector of the cell a. In this framework, the acoustic sum rule writes [START_REF] Esfarjani | Method to extract anharmonic force constants from first principles calculations[END_REF][START_REF] Esfarjani | Heat transport in silicon from first-principles calculations[END_REF] :

i (1) Θ α i (0) = 0 ∀ (α) jb (2)
Θ αβ ij (0, b) = 0 ∀ (i, αβ) kc (3)
Θ αβγ ijk (0, b, c) = 0 ∀ (ijb, αβγ) ld (4)
Θ αβγδ ijkl (0, b, c, d) = 0 ∀ (ijkbc, αβγδ) with (p) Θ α 1 ...αp i 1 .
..ip (a 1 , ..., a p ) the effective IFC at the p th order corresponding to the interaction between the (i th 1 , ..., i th p ) atoms in cells (a 1 , ..., a p ), respectively. These equations lead to some relations between all the coefficients of each IFC. By using some mathematics (and contracting the summation index over atoms), it can be showed that Equation 3 can be rewritten as :

F α i = - (1)
Θ α i - β,j =i (2)
Θ αβ ij (u β j -u β i ) - 1 2 βγ j =i,k =i (3) Θ αβγ ijk (u β j -u β i )(u γ k -u γ i ) - 1 6 βγδ j =i,k =i,l =i (4) Θ αβγδ ijkl (u β j -u β i )(u γ k -u γ i )(u δ j -u δ i ) (9)
We can easily see that, at each order, the on-site coefficients are not involved in the calculation of the total force and that the invariance of the system under arbitrary translation is fully achieved. The constraints coming from the acoustic sum rules lead to a reduction of the IFC coefficients and we will show in Section 2.7 how to proceed.

The rotational invariance

The system (force and energy) has to remain invariant after a rotation. This imposes some relations between the IFC of different orders [START_REF] Esfarjani | Method to extract anharmonic force constants from first principles calculations[END_REF][START_REF] Esfarjani | Heat transport in silicon from first-principles calculations[END_REF] :

0 = i,αβ (1)
Θ α i (0)τ β i αβν ∀ (ν) 0 = jb,βγ (2)
Θ αβ ij (0, b)τ γ j βγν + β (1)
Θ β i (0) βαν ∀ (i, αν) 0 = γ (2) Θ γβ ij (0, b) γαν + γ (2)
Θ αγ ij (0, b) γβν + kc,γδ (3)
Θ αβγ ijk (0, b, c)τ δ k γδν ∀ (ijb, αβν) 0 = δ (3) Θ δβγ ijk (0, b, c) δαν + δ (3) Θ αδγ ijk (0, b, c) δβν + δ (3) Θ αβδ ijk (0, b, c) δγν + ld,δµ (4)
Θ αβγδ ijkl (0, b, c, d)τ µ l δµν ∀ (ijkbc, αβγν)
with αβγ the Levi-Civita symbols. As for the acoustic sum rules defined in the previous subsection, one explains in Section 2.7 how to impose such constraints.

The crystal symmetries (I)

Using the crystal symmetry operations S, it is possible to pursue the reduction process of the IFC coefficients:

• by deducing all the IFC p connected to a reference IFC ref p through a crystal symmetry (that is the subject of this section). In this case, the symmetries S used are the space group operations (mirrors, rotations, translations including non-symmorphic ones). • by determining the nonzero coefficients really independent in the IFC ref p (this will be the subject of the next section). In this case, the symmetries S used are the point group operations.

By means of space group operations, it is possible to deduce all the interactions from a set of reference ones. Schematically, this reduction is achieved using simultaneously two kinds of operations : by taking into account the translations, it reduces the complexity of the supercell to the first unitcell at most, and by considering the mirrors and rotations, it is possible to build some shells of interaction (see Figure 1) around reference atoms in the first unitcell. Therefore, using the symmetry operation S, each "image" interaction, with the atoms indexed by the S(i)...S(l) letters, can be deduced from a "reference" interaction involving the atoms i...l, such as :

(1) Θ α S(i) = µ S αµ (1)
Θ µ i ( 10
)
(2)

Θ αβ S(i)S(j) = µν S αµ S βν (2)
Θ µν ij (11)
(3)

Θ αβγ S(i)S(j)S(k) = µνξ S αµ S βν S γξ (3)
Θ µνξ ijk ( 12
) (4) Θ αβγδ S(i)S(j)S(k)S(l) = µνξo S αµ S βν S γξ S δo (4)
Θ µνξo ijkl ( 13
)
The number of coefficients falls down significantly after this work. If 10 shells of pair interactions are present, there are at worst 3 2 *10=90 coefficients to determine for the second order... and no longer (3*N a ) 2 =90 000 as before.

The crystal symmetries (II)

The last symmetry rule is probably the most subtle and tedious to carry out. The aim of this subsection is to obtain a priori the exact number of nonzero coefficients really independent in the reference matrices IFC ref p . The symmetry operations S used here belong to the point group, which forms a subgroup of the orthogonal group. Consequently, all the symmetry operations S (in the cartesian space) belonging to this group have the following properties :

• The inverse equals its transpose : S -1 = S T

• The norm of the (three) eigenvalues equals to 1: λ = (±1, ±1, ±1) when S is symmetric, or λ = (±1, e iθ , e -iθ ) when S is antisymmetric.

• The norm of the (three) eigenvectors equals to 1: |p| = 1 These properties will be used thereafter and will strongly simplify the development. We have to consider two peculiar cases : some symmetry operations S may let the interaction invariant and/or some others may reverse the interaction. This imposes some constraints on the IFC and thus reduces the number of coefficients. Note that the symmetry operations leaving the interaction invariant (or reversed) form a subgroup of the point group, so a large set of constraints may consistently act on the IFC. In the following, we restrict the development to the second order. The reasoning is the same for all the others.

If the symmetry leaves the bond invariant.... the IFC ref 2

is not affected by the transformation, so S(i) = i and S(j) = j in Equation 11 which becomes :

(2)

Θ αβ ij = µν S αµ S βν (2)
Θ µν ij
We can apply the eigenvectors of the S symmetry operation and obtain :

αβ p α r p β s (2)
Θ αβ ij = λ r λ s µν p µ r p ν s (2) Θ µν ij ∀(r, s = 1, 2, 3)
By gathering the left and right members we get:

(λ r λ s -1) αβ p α r p β s (2)
Θ αβ ij = 0 ∀(r, s = 1, 2 , 3)
The left hand-side of the equation equals to zero if i) the term between brackets is null or if ii) the summation is null. The eigenvalues of S being equal to ±1 or e ±iθ , it happens very frequently that their product equals to 1. Consequently, the condition i) and the previous equality is naturally fulfilled in numerous cases. Otherwise, we have :

αβ p α r p β s (2)
Θ αβ ij = 0 if λ r λ s = 1 ∀(r, s = 1, 2, 3) (14)
These equations impose some constraints on the IFC ref

2

(linear relation between the coefficients). We will see in the following how to extract the nonzero independent coefficients.

If the symmetry reverses the bond.... the permutation also affects IFC ref 2 , so S(i) = j and S(j) = i in Equation 11 which becomes :

(2)

Θ αβ ji = µν S αµ S βν (2)
Θ µν ij Using the symmetry properties of the IFC then by applying the eigenvectors and gathering the left and right hand-side, one obtains :

αβ (λ r λ s p α r p β s -p β r p α s ) (2)
Θ αβ ij = 0 ∀(r, s = 1, 2, 3) (15)
As already highlighted when the interaction is kept invariant, these equations also impose some constraints if the term within brackets is not null.

All the symmetries together

For each order and shell independently, all the equations (see equations 14 and 15) are collected then orthonormalized, resulting in a system of Λ p independent homogeneous linear equations (at the p th order). This system may be rewritten as :

3 p κ=1 Υ κ λp (p) Θ κ i 1 ...ip = 0 with λ p = 1, ...Λ p
with the superscripts α 1 ...α p contracted in κ, and the matrix Υ κ λp coming from all the nonzero terms between brackets. The solutions of such system is obtained by computing the kernel of the application. The result has the following form : By taking into account all the interactions of each shell at each order, we finally recover Equation 5:

F α i (t) = pλp f α i,pλp (u(t))θ pλp .
All the terms depending on the displacements u α i (t) and symmetries S αβ are hidden in the f matrix. We obtain a system of linear equations with the nonzero coefficients of the effective IFC {θ pλp } as the unknown variables. Its expanded form is :

           F x 1 (1) F y 1 (1) F z 1 (1)
. . .

F z Na (1)
. . .

F z Na (N t )            =             f x 1,11 (1) f x 1,12 (1) . . . f x 1,pΛp (1) f y 1,11 (1) f y 1,12 (1) . . . f y 1,pΛp (1) f z 1,11 (1) f z 1,12 (1) . . . f z 1,pΛp (1)
(N t ) f z Na,12 (N t ) . . . f z Na,pΛp (N t )                        θ 11 θ 12 θ 13 . . . θ 1Λ 1 . . . θ pΛp           
We stress that the f matrix is really huge. At the second order, its size is

(3 * N a * N t ) × Λ 2 .
If we consider N t =1000 time steps, N a =100 atoms and Λ 2 =20 coefficients, this matrix has 300 000 lines and 20 columns.

The constrained least squares problem

The whole set of coefficients {θ pλp } can be obtained by reversing the previous equation. An elegant way [START_REF] Esfarjani | Method to extract anharmonic force constants from first principles calculations[END_REF][START_REF] Hellman | Lattice dynamics of anharmonic solids from first principles[END_REF] to perform such reversion is to compute the Moore-Penrose pseudoinverse (marked with a † superscript in the following, see Equation 8):

   θ 11 . . . θ pΛp    =    f x 1,11 (1) . . . f x 1,pΛp (1)
. . . . . . . . .

f z Na,11 (N t ) . . . f z Na,pΛp (N t )    † .    F x 1 (1)
. . .

F z Na (N t )   
The usual and mathematical definition of a pseudoinverse Γ † is :

Γ † = Γ T (Γ.Γ T ) -1 (16)
This definition implies to carry out a reversion of the f matrix. In practice, this method is proved to be numerically unstable and the computational cost turn out to be really expensive, due to the size of the matrix. A more simple and accurate way is possible, by performing a singular value decomposition :

Γ † = V.Σ † .U T (17)
with V and U some square matrices and Σ a diagonal matrix. Some linear algebra routines (dgesdd) enable to perform such decomposition easily. The Σ being diagonal, this one can be conveniently reversed, which makes this method much more stable and less expensive than the previous one.

We have seen previously that the acoustic sum rules and the rotational invariance of the system impose some constraints on the IFC. In the same way, we will see in the following that the dynamical matrix and the elastic tensor also impose other constraints on IFC. If one rewrites these Q linear constraints, function of the θ i ≡ θ pλp coefficients, in a matrix C, we obtain an homogeneous linear equation system :

pΛp i=1 C T q,i θ i = 0 with q = 1, ..., Q (18)
In order to satisfy these constraints we have to :

1. minimize S = ||F M D -f .Θ|| 2 , as previously performed by least squares method with Eq.7

2. subject to C.Θ = 0.

That can be done by building the Lagrangian function L(θ, ζ) with ζ the Lagrange multipliers :

L(θ, ζ) = S(θ) + Q q=1 ζ q pΛp i=1 C T q,i θ i (19)
The optimality conditions are :

∂L ∂θ i = 2 pΛp j=1 (f T .f ) ij θ j -2(f T .F M D ) i + Q q=1 ζ q C q,i = 0 ∂L ∂ζq = Q q=1 C T q,i θ i = 0
In matrix-vector form, these equations are more readable and are written :

2f T .f C T C 0 . Θ ζ = 2f T .F M D 0 ( 20
)
As already carried out for the least square problem, the constrained least squares (CLS) solution can be obtained by computing the pseudoinverse of the first matrix in the left hand-side member. Note that the matrix to invert is smaller than the one shown previously. At the second order, its size is

(Λ 2 + Q) × (Λ 2 + Q).
This one no longer depends on the number of time steps and atoms. For Λ 2 =20 coefficients and Q=5 constraints this matrix has only 25 lines and columns.

In conclusion, we have shown in this section how to strongly reduce the number of IFC coefficients using symmetries and invariances of the system. In this implementation, the CLS method used to compute the IFC coefficients is applied only once (for all orders simultaneously) and all the IFC coefficients are obtained together. Consequently, the 2 nd order effective IFC may change depending on whether the 3 rd and/or 4 th order is included or not in the CLS resolution. This strategy could be problematic if the Taylor expansion is divergent (when the 3 rd and/or 4 th order become larger than the 2 nd order). One way to overcome this break down is to solve each order, one after eachother. The procedure is to compute the residual of the forces at the n th order, then apply the CLS method and evaluate the (n + 1) th effective IFC. In this case the 2 nd order effective IFC no longer changes depending on whether highest order are included or not in the expansion. In a future version of this code, we will introduce this new strategy and let the user choose to solve all the orders simultaneously or successively.

What are the effective IFC useful for?

In this section we show that the knowledge of effective IFC allows us to calculate a large number of physical quantities that characterize a material: phonon spectrum, free energy, elastic constants and moduli, sound speed, Grüneisen parameter, thermal expansion... Since the effective IFC are temperature dependent, all the quantities resulting from them exhibit an explicit dependence with respect to temperature.

The phonon spectrum

In order to evaluate the phonon spectrum and the vibrational density of states, we need to split up again the index over atoms as perform in the subsection 2.2. Once the

(2)
Θ αβ ij (a, b) tensor is built (after the Moore-Penrose), it is quite direct to compute the phonons modes of the system by performing its Fourier Transform. We obtain the so-called dynamical matrix :

D αβ ij (q) = b (2) Θ αβ ij (0, b) M i M j exp (iq.[R(b) -R(0)]) (21)
with q a wave vector belonging to the Brillouin zone (BZ) and M i the mass of atom i. Such equations of motion can be diagonalized in order to obtain the phonon eigen frequencies ω s (q) and eigen modes X α is (q) of the lattice :

β,j D αβ ij (q)X β js (q) = ω 2 s (q)X α is (q) ( 22
)
In this equation (representative of an harmonic oscillator) the quantum number s refers to the quantification of the eigen frequencies, with 3 * N a energy levels at all. Note that the phonon spectrum depends on the presence/absence of the higher terms in the model.

The vibrational Density Of States (vDOS) g(ω) is defined by :

g(ω) = 1 3N a 3Na s=1 q∈BZ δ(ω -ω s (q)) such as ωmax 0 g(ω)dω = 1 (23)
with ω max the highest phonon frequency of the system.

Thermodynamic quantities

When one wants to evaluate the free energy F (V, T ) of a system as a function of the temperature and the volume, the ususal process is to split up the cold and vibrational contributions such as :

F (V, T ) = U 0 (V ) + F vib (V, T ) (24)
with U 0 (V ) the zero temperature energy of the system and F vib (V, T ) the vibrational free energy. If the cold contribution can be estimated through an independent ground state calculation, this one can also be computed as the average over time steps of the difference between the the total energy U MD (t) coming from the AIMD simulation and the expansion of the energy as defined by Equation 2:

U 0 (V ) = U MD (t) - P p=1 1 p ! α 1 ...αp i 1 ...ip (p) Θ α 1 ...αp i 1 ...ip p k=1 u α k M D,i k (t) (25) with = 1
Nt Nt t=1 . If there is no strong modification of the equilibrium structure of the crystal due to temperature effects, these two approaches have to give similar results.

Concerning the vibrational part of the free energy F vib (see Equation 24) but also in a similar way for the internal energy U vib , the constant-volume specific heat C vib,V = ∂U ∂T V and the entropy S vib = -∂F vib ∂T V , they can be evaluated in the framework of the three-dimensional quantum harmonic crystal using the vDOS [START_REF] Lee | Ab initio calculation of the thermodynamic properties and atomic temperature factors of sio 2 α-quartz and stishovite[END_REF] :

F vib = 3N a β ωmax 0 ln 2 sinh( β ω 2 ) g(ω)dω (26)
U vib = 3N a ωmax 0 ω 2 coth( β ω 2 )g(ω)dω (27)
C vib,V = 3N a k B ωmax 0 β ω 2 sinh( β ω 2 ) 2 g(ω)dω (28)
S vib = 3N a k B ωmax 0 β ω 2 coth( β ω 2 ) -ln 2 sinh( β ω 2 ) g(ω)dω (29)
with β = k B T . Note that these formulas (and the others thereafter in this section) are the common harmonic expressions and have to be modified in order to include the corrections coming from the third and fourth orders [START_REF] Leibfried | Theory of anharmonic Effects in Crystals[END_REF][START_REF] Barron | Thermal and elastic properties of crystals at low temperatures[END_REF][START_REF] Cowley | Zero sound, first sound and second sound of solids[END_REF][START_REF] Cowley | Anharmonics crystals[END_REF][START_REF] Shukla | Atomic mean-square displacement of a solid : A green's-function approach[END_REF][START_REF] Allen | Anharmonic phonon quasiparticle theory of zeropoint and thermal shifts in insulators: Heat capacity, bulk modulus, and thermal expansion[END_REF][START_REF] Oba | First-principles study of phonon anharmonicity and negative thermal expansion in scf 3[END_REF]. This work is in progress and will be achieved for a next release of the code. In a similar way [START_REF] Minakov | Melting curves of metals with excited electrons in the quasiharmonic approximation[END_REF], one can also compute the mean square displacement :

< u 2 >= 3 2M a ωmax 0 coth( β ω 2 ) g(ω) ω dω (30)
All these thermodynamic quantities display three kinds of dependence with respect to the temperature. The first one is given by β and originates from the Bose-Einstein statistic and corresponds to the filling of the energy levels (including the zero-point energy), another one (ω(V (T ))) more implicit is provided by the thermal expansion of the volume, and a last one (ω(T )) is due to the explicit dependence of the phonon frequencies with respect to the temperature. If the first one is taken into account in a quantum harmonic approach, the second one appears when considering the QHA approximation, and the third one comes from the explicit treatment of the temperature carried out in this work.

Elastic constants and moduli

The effective IFC enable to obtain almost all the thermal elastic properties of a material. The elastic constants are defined using the Hooke law which relies on the strain tensor αβ to the stress one σ γδ , through the elastic tensor C αβγδ :

σ αβ = C αβγδ γδ (31)
The most famous way employed to compute the elastic tensor is the finite differences method. This approach, coupled with AIMD simulations, enable to take into account temperature effects. But numerous simulations have to be performed and the computational time can be prohibitive (in order to have long AIMD trajectories). Another approach is to compute phonons as a function of temperature and to evaluate the elastic constants using the slopes of the acoustic branches. This one is connected to a more elegant formulation using the IFC [START_REF] Leibfried | Theory of anharmonic Effects in Crystals[END_REF][START_REF] Wallace | Statistical physics of crystals and liquids[END_REF]. Following D.C. Wallace [START_REF] Wallace | Statistical physics of crystals and liquids[END_REF], let us define the A αβγδ quantity :

A αβγδ = 1 2V ij (2)
Θ αβ ij d γ ij d δ ij with d γ ij = τ γ i -τ γ j (32)
The elastic constants are then obtained by using this simple following equation :

C αβγδ = A αγβδ + A βγαδ -A αβγδ (33)
Using the Voigt notation and the Voigt formulation, we can define the isothermal Bulk K T and Shear G moduli (but also the Young modulus E and the Poisson's ratio ν) as follows : [START_REF] Paulatto | First-principles calculations of phonon frequencies, lifetimes, and spectral functions from weak to strong anharmonicity: The example of palladium hydrides[END_REF] It is important to note that this method needs to have long range IFC in order to have converged quantities [START_REF] Schnell | Tight-binding calculations of the elastic constants and phonons of hcp zr: Complications due to anisotropic stress and long-range forces[END_REF].

K T = ((C 11 + C 22 + C 33 ) + 2(C 12 + C 13 + C 23 ))/9 (34) G = ((C 11 + C 22 + C 33 ) -(C 12 + C 13 + C 23 ) + 3(C 44 + C 55 + C 66 ))/15

Grüneisen parameter & thermal expansion

As well known, the thermal expansion α V = 1 V ∂V ∂T P , the thermodynamic Grüneisen parameter γ = V ∂P ∂U V , the isothermal compressibility 34) and the specific heat C V (see Equation 28) are connected to each other through the following equation:

β T = -1 V ∂V ∂P T (= 1/K T , see Equation
α V = γC V β T V (36)
In order to compute γ in the lattice dynamic framework, we introduce the mode Grüneisen parameter γ s (q) and the specific heat C V,s (q) per phonon mode ω s (q):

γ s (q) = - V ω s (q) ∂ω s (q) ∂V T ( 37
)
C V,s (q) = k B β ω s (q) 2 sinh( β ωs(q) 2 ) 2 (38)
It can be shown that the thermodynamic Grüneisen parameter γ is written :

γ = s,q γ s (q)C V,s (q) C V (39)
If the system undergoes an anisotropic strain γδ , such as κ 2 κ = 1 in the Voigt notation, the displacement of atom k becomes u γ k = δ γδ τ δ k . In this framework, the mode Grüneisen parameter rewrites [START_REF] Lee | Anharmonic interatomic force constants and thermal conductivity from grüneisen parameters : An application to graphene[END_REF] :

γ γδ s (q) = - 1 3ω s (q) ∂ω s (q) ∂ γδ T = - 1 6ω 2 s (q) ∂ω 2 s (q) ∂ γδ T (40)
Using Equations 22 and 21, we obtain :

γ γδ s (q) = - 1 6ω 2 s (q) ∂ ∂ γδ ij,αβ X α is (q) b 1 M i M j ∂ 2 U ∂u α i ∂u β j exp (iq.R(b)) X β js (q)
(41) In this equation one has to compute the third derivative with respect to displacements and strain. This one can be evaluated using the expansion of the Hamiltonian up to the third order:

∂ ∂ γδ ∂ 2 U ∂u α i ∂u β j = ∂ ∂ γδ (2) Θ αβ ij + kγ (3) Θ αβγ ijk u γ k = k (3) Θ αβγ ijk τ δ k ( 42
)
We obtain the common definition of the mode Grüneisen parameter as a function of the interatomic force constants.

γ γδ s (q) = - 1 6ω 2 s (q) ijk,bc,αβ (3)
Θ αβγ ijk (0, b, c)

X α is (q)X β js (q) M i M j τ δ k exp [iq.R(b)] (43)
Once these parameters obtained, the thermodynamic Grüneisen parameter γ (see Equation 39) and the thermal expansion α V (see Equation 36) can be computed directly. In the case of an anisotropic crystal, the thermal expansion matrix rewrites :

α αβ V = γδ S αβγδ qs C V,s (q)γ γδ s (q) ( 44
)
with S αβγδ ≡ C -1 αβγδ the compliance matrix (see Equation 33). The volumetric thermal expansion α V can be obtained by summing up the linear thermal expansion coefficients α αα V arranged along the diagonal of the thermal expansion matrix : α V = i α αα V .

3.5. Thermal pressure, isentropic quantities, sound speed & Debye temperature If the electronic entropy contribution is not considered, the total stress tensor can be split into two parts (see Equation 24) :

σ αβ = 1 V ∂F ∂ αβ T = σ 0,αβ + σ vib,αβ (45)
with the first part coming from the cold energy U 0 (V ) and the second one being the thermal stress contribution deriving from the thermal free energy (see Equation 26) as :

σ vib,αβ = - s,q γ αβ s (q) U vib,s (q) V with ( 46
)
U vib,s (q) = ω s (q) 2 coth( β ω s (q) 2 ) ( 47
)
This result is really important since it becomes possible to verify that the pressure (P = -1 3 Tr[σ]) computed using the AIMD code is well reproduced by the one computed using a-TDEP. The latter part involving the mode Grüneisen parameter, it is a good way to check its convergence.

The thermodynamic Grüneisen parameter and the thermal expansion are also useful to compute the isentropic (or adiabatic) compressibility β S and the constant-pressure heat capacity C P :

β S = β T (1 + α V γT ) = 1 K S ( 48
)
C P = C V (1 + α V γT ) ( 49
)
Note that β T and C V are already known (see Equations 34 and 28). Using the isentropic compressibility, the longitudinal v l and transverse v t sound speeds can be also obtained :

v l = K S + 4G/3 ρ and v t = G ρ (50)
and the Debye temperature Θ D becomes [START_REF] Shao | First-principles study on the elastic properties of cu2gese3[END_REF][START_REF] Ashcroft | Solid State Physics[END_REF] :

Θ D = k B 6π 2 N a V 1/3 v m with v m = 1 3 2 v 3 t + 1 v 3 l -1/3 (51)

a-TDEP inside the Abinit package

The algorithms and all the thermodynamic quantities detailed in the previous sections have been implemented in the Abinit package. If the a-TDEP code is present in Abinit since the Version 8.6 (released on November 3, 2017), the current described version is implemented in the Version 8.12.

As deeply integrated with Abinit, it follows all its coding rules. It uses Fortran2003 language and makes use of data structures close to a Object-Oriented Programming. Thanks to the structures and features already implemented in Abinit, many functionalities are directly taken from it for easiness, reliability and robustness. Regarding parallelism, even if Abinit has several levels of parallelism [START_REF] Bottin | Large-scale ab initio calculations based on three levels of parallelization[END_REF] with both the Message Passing Interface (MPI) and OpenMP (threads), a-TDEP does not use parallelism yet. Nevertheless, compiling with a multithreaded linear algebra library will increase the performances of the code.

As a-TDEP takes place in the source tree of Abinit, the build system of it is fully compatible and nothing more than the build system of Abinit is needed. Actually, during the configure step and build step of Abinit, everything is set for a-TDEP. Therefore all the options and tuning of Abinit can be used inside a-TDEP like tuning the linear algebra library and so on without any effort [START_REF] Gonze | Abinit: First-principles approach to material and nanosystem properties[END_REF]. At the end of the build, a new executable is available (as Abinit and many other post-process) and can be directly used.

For launching the executable, a-TDEP follows the same strategy as Abinit. A files file can be provided (or a prompt will be displayed) containing :

1. The input file (input.in) 2. The molecular dynamic trajectory (the HIST.nc file provided by Abinit in the NetCDF format)

3. A root base name for all the output files (root).

Many output files are generated by the code, mainly for checking purpose.

The main output file is :

1.

[root].out. It includes an echo of the input variables, some intermediary results, the definition of the various shells of interaction, the IFC for all the atoms in each shell, the elastic constants and moduli, the energy of the model...

Other important files are also obtained :

2.

[root]omega.dat contains the dispersion of phonon frequencies (in meV) along a path in the Brillouin Zone,

3.

[root]thermo.dat lists all the thermodynamic quantities obtained by considering the system as a quantum harmonic crystal : specific heat, vibrational energy, entropy and free energy. It also gives all these contributions as a function of temperature in the harmonic approximation.

4.

[root]phonons-band.yaml containing the phonon band structure in the Phonopy convenient format [START_REF] Togo | First-principles calculations of the ferroelastic transition between rutile-type and cacl 2 -type sio 2 at high pressure[END_REF] 5.

[root]sym.dat details all the symmetry operations of the bravais lattice,

6.

[root]qpt.dat defines the q-point grid used to compute the phonon frequencies contained in the omega.dat file,

7.

[root]xred average.xyz includes the ideal and average positions in the supercell, 8.

[root]Indym*.dat contain all the symmetry relations between one or two atoms in the unitcell or the supercell, 9.

[root]vdos.dat displays the vibrational density of states (in meV),

10.

[root]dij.dat lists the dynamical matrices for a particular set of qpoints,

11.

[root]etotMDvsTDEP [p].dat compares the MD energies with the one computed using the p th order effective IFC as a function of time step (these ones must be superimposed, as much as possible) ,

12.

[root]fcartMDvsTDEP [p].dat plots the MD forces wrt the forces computed using the p th order IFC (the cloud of points must be closer to the first bisector),

13.

[root]eigenvectors.dat lists all the eigenvectors for a particular set of q-points,

14.

[root]nbcoeff-phij.dat/[root]nbcoeff-psij.dat show how the number of IFC coefficients are reduced (for each shell and each symmetry) at the second and third order, respectively.

15. ... For reliability and robustness, a-TDEP benefits from the full Abinit test farm. Abinit already includes more than 1200 tests that enables not only a strong guarantee on the results but also on portability of the code. a-TDEP takes advantages of all this pre-existent framework and includes its owns tests on several materials, geometry and physical quantities, including those presented in the following. This makes the code safer and more reliable.

To help the user, an external open source tool is also available. This one is called qTdep and is interfaced with the present version of a-TDEP. Its use is very simple. First provide the set of configurations or molecular dynamic trajectory (first tab of Figure 2). After that, qTdep will try to figure out what the unit cell is and if so, it will compute the multiplicity. If it can not determine it, the user has to provide an input file of this reference structure. Every time the unitcell is modified, the multiplicity of the supercell is recalculated. Finally, the last tabulation on Figure 2 shows the options available to perform the calculation. The order of the calculation (and thus the properties that can be extracted) the cutoff radius for each order and some boolean options to tweak the calculation. Then, when the changes are applied, the calculation is launched (a-TDEP has to be in the PATH environment variable) and when the calculation succeeds, the phonon dispersion curves are plotted. The input file for a-TDEP is created automatically by qTdep so the user can change it by hand and tweak even more the calculation.

Examples illustrating the a-TDEP capabilities

In this section, we are interested in three materials : Si, Zr and U. The results provided below are not new or original. This section is just the occasion (i) to put into practice a large number of the a-TDEP capabilities (phonon frequencies, thermal expansion, Grüneisen parameters, elastic constants...) and (ii) to demonstrate that our home made implementation reproduce some results already published by Hellman and coworkers. )

R cut3 =8 a.u. R cut3 =10 a.u. R cut3 =12 a.u. R cut3 =14 a.u. R cut3 =16 a.u.

Negative thermal expansion in silicon

The negative thermal expansion (NTE) of Si between 20 and 120 K (see Fig. 3) is a longstanding issue in condensed matter physics [START_REF] Batchelder | Lattice constants and thermal expansivities of silicon and of calcium fluoride between 6 deg and 322 deg k[END_REF][START_REF] Carr | Thermal expansion of germanium and silicon at low temperatures[END_REF][START_REF] Ibach | Thermal expansion of silicon and zinc oxide (i)[END_REF][START_REF] Nilsson | Study of the homology between silicon and germanium by thermal-neutron spectrometry[END_REF][START_REF] Slack | Thermal expansion of some diamondlike crystals[END_REF][START_REF] Weinstein | Raman scattering and phonon dispersion in si and gap at very high pressure[END_REF][START_REF] Lyon | Linear thermal expansion measurements on silicon from 6 to 340 k[END_REF][START_REF] Okada | Precise determination of lattice parameter and thermal expansion coefficient of silicon between 300 and 1500 k[END_REF][START_REF] Strauch | Phonon eigenvectors in si determined by inelastic neutron scattering[END_REF][START_REF] Kulda | Inelastic-neutronscattering study of phonon eigenvectors and frequencies in si[END_REF][START_REF] Rignanese | Ab initio study of the volume dependence of dynamical and thermodynamical propoerties of silicon[END_REF][START_REF] Fabian | Thermal expansion and grüneisen parameters of amorphous silicon: A realistic model calculation[END_REF][START_REF] Broido | Lattice thermal conductivity of silicon from thermal interatomic potentials[END_REF][START_REF] Esfarjani | Method to extract anharmonic force constants from first principles calculations[END_REF][START_REF] Esfarjani | Heat transport in silicon from first-principles calculations[END_REF][START_REF] Härkonen | Ab initio lattice dynamical studies of silicon clathrate frameworks and their negative thermal expansion[END_REF][START_REF] Tadano | Anharmonic force constants extracted from first-principles molecular dynamics: applications to heat transfer simulations[END_REF][START_REF] Kim | Phonon anharmonicity in silicon from 100 to 1500 k[END_REF][START_REF] Middelmann | Thermal expansion coefficient of single-crystal silicon from 7 k to 293 k[END_REF][START_REF] Argaman | Ab initio study of the phononic origin of negative thermal expansion[END_REF]. If a good description of such unusual property can be provided in the framework of QHA [START_REF] Rignanese | Ab initio study of the volume dependence of dynamical and thermodynamical propoerties of silicon[END_REF][START_REF] Argaman | Ab initio study of the phononic origin of negative thermal expansion[END_REF][START_REF] Kim | Nuclear quantum effect with pure anharmonicity and the anomalous thermal expansion of silicon[END_REF], we will show in the following that a quantitative agreement is only recovered when thermal effects are explicitly taken into account, as shown by Kim et al. [START_REF] Kim | Nuclear quantum effect with pure anharmonicity and the anomalous thermal expansion of silicon[END_REF].

In order to compute accurate effective IFC by means of the least squares method, we have to accumulate a sufficiently large set of forces and positions. So, we perform AIMD simulations during 10 ps and select 250 configurations within (in the thermalized part of the trajectory). The supercell contains 216 atoms and the calculations are carried out in the NVT ensemble, at room temperature and using the experimental lattice parameter (a exp =5.43 Å) rather than the theoretical one (a th =5.40 Å).

Our calculations are performed by using the computational code ABINIT [91] over hundreds to thousands of processors [START_REF] Bottin | Large-scale ab initio calculations based on three levels of parallelization[END_REF]. A norm-conserving pseudopotential is generated following the Troullier-Martins scheme [START_REF] Troullier | Efficient pseudopotentials for plane-wave calculations[END_REF] in the Kleinman-Bylander form [START_REF] Kleinman | Efficacious form for model pseudopotentials[END_REF], considering 4 electrons in the valence and a cutoff radius of 2.1 bohr leading to a cutoff energy equal to 380 eV. Exchange and correlation energy was treated with the local density approximation LDA using the Teter parametrization [START_REF] Goedecker | Separable dual-space gaussian pseudopotentials[END_REF]. A careful treatment of the electronic density integration has been achieved with the use of a (2×2×2) Monkhorst-Pack mesh [START_REF] Monkhorst | Special points for brillouin-zone integrations[END_REF].

Two important parameters drive the number of IFC coefficients really computed during the a-TDEP calculation : R cut2 and R cut3 , the cutoff radii at the 2 nd or 3 rd order, respectively. We are tempted to put the highest possible value for these parameters. However, we face two limitations. The first one is the size of the supercell; the cutoff radius cannot exceed half of its lowest dimension. The second one is the computational cost : the larger these parameters, the greater the number of IFC coefficients to calculate. Therefore, five values has been used to check the convergence : 8, 10, 12, 14 and 16 a.u.. These cutoff radii lead to treat 6, 10, 16, 27 and 37 then 36, 95, 219, 711 and 1408 independent IFC coefficients, at the 2 nd order and 3 rd order respectively.

In Fig. 3, we display the variation of the linear thermal expansion α with respect to the temperature. In this figure, we compare the results obtained for the five values of the R cut3 parameter (with R cut2 =16 a.u.) with five sets of experimental data [START_REF] Carr | Thermal expansion of germanium and silicon at low temperatures[END_REF][START_REF] Ibach | Thermal expansion of silicon and zinc oxide (i)[END_REF][START_REF] Slack | Thermal expansion of some diamondlike crystals[END_REF][START_REF] Lyon | Linear thermal expansion measurements on silicon from 6 to 340 k[END_REF][START_REF] Middelmann | Thermal expansion coefficient of single-crystal silicon from 7 k to 293 k[END_REF] (here, we do not consider the results of Batchelder and Simmons [START_REF] Batchelder | Lattice constants and thermal expansivities of silicon and of calcium fluoride between 6 deg and 322 deg k[END_REF] since we find that there are inconsitent with the other experiments). As a function of the 3 rd order cutoff, we obtain two groups of results : a poor agreement with experiments when the R cut3 parameter is "too low" (for 8 and 10 a.u.) and a better one when the R cut3 parameter is larger (for 12, 14 and 16 a.u.). At R cut3 =16 a.u., the agreement is surprisingly perfect, especially around the temperature used in the AIMD simulations (300 K), even if the convergence seems not completely achieved as a function of the cutoff radius. Our results are also in excellent agreement with the ones obtained using the original implementation of O. Hellman [START_REF] Kim | Nuclear quantum effect with pure anharmonicity and the anomalous thermal expansion of silicon[END_REF]. At higher temperature, the inclusion of the 4 th order terms in the expansion of the energy (not considered here) are needed [START_REF] Yue | Metric for strong intrinsic fourth-order phonon anharmonicity[END_REF] to conserve a good agreement with experiments [START_REF] Kim | Nuclear quantum effect with pure anharmonicity and the anomalous thermal expansion of silicon[END_REF].

The convergence of the thermal expansion is in fact driven by the three quantities involved in its calculation : the thermodynamic Grüneisen parameter, the Bulk Modulus and the specific heat (see Equation 36). In Figure 4, we show the variation of the elastic constants C 11 , C 12 and the Bulk Modulus K T as a function of the 2 nd Order cutoff radius R cut2 (R cut3 =0 in this case). We see that the convergence of these three moduli is very long as a function of the cutoff radius, as already exposed for hcp Zr [START_REF] Schnell | Tight-binding calculations of the elastic constants and phonons of hcp zr: Complications due to anisotropic stress and long-range forces[END_REF] or TiN [START_REF] Steneteg | Temperature dependence of tin elastic constants from ab initio molecular dynamics simulations[END_REF]. Concerning the Bulk Modulus, the agreement with experiments is very good for R cut2 =8 or 10 a.u., due to an error compensation (negative error for C 11 and positive error for C 12 , see Eq. 34), then very bad for R cut2 =12 or 14 a.u. due to a better description of the C 12 , and finally gives an intermediary value, with 7% of error wrt the experimental one.

Another important contribution to the thermal expansion is the thermodynamic Grüneisen parameter γ, built using the mode Grüneisen parameter γ s (q) (see Eq. 39). In Fig. 5, we can see a small variation in their dispersion curves when moving from the R cut3 =8 or 10 a.u. first group to the R cut3 =12, cutoff radius is significant (in particular for the acoustic branches). Once again, two groups can be distinguished, according to the inclusion (or not) of the {100} shell in the a-TDEP calculation : a first group represented by the R cut2 =8 a.u. radius and a second one represented by the R cut2 =12 a.u. radius. For R cut2 =R cut3 =16 a.u., the agreement with experiments [START_REF] Nilsson | Study of the homology between silicon and germanium by thermal-neutron spectrometry[END_REF][START_REF] Kulda | Inelastic-neutronscattering study of phonon eigenvectors and frequencies in si[END_REF] is excellent. Only the optical branches are slightly lower than experimental data (as also obtained by other authors before [START_REF] Tadano | Anharmonic force constants extracted from first-principles molecular dynamics: applications to heat transfer simulations[END_REF][START_REF] Esfarjani | Method to extract anharmonic force constants from first principles calculations[END_REF]). The use of the theoretical lattice parameter rather than the experimental one slightly increases these optical branches, but it is not sufficient to have a perfect agreement.

Γ X Γ L -3 -2 -1 0 1 2 Mode Gruneisen parameter R cut3 =8 a.u. R cut3 =12 a.u. Exp
(THz) R cut2 =8 a.u. (V exp ) R cut2 =12 a.u. (V exp ) R cut2 =R cut3 =16 a.u. (V exp ) R cut2 =R cut3 =16 a.u. (V th ) Exp1 Exp2
At last, we present some very useful quantities which are available as outputs of a-TDEP and allow to check the fit procedure. After applying the least square method, the model (TDEP) has to reproduce correctly the AIMD simulations. In Fig. 7, we plot the TDEP forces F T DEP as a function of the AIMD ones F AIM D and, in the two insets, the energies and their differences (for the 250 configurations selected at the beginning). We see clearly that the inclusion of the 3 rd Order improves the fit. The standard deviation on forces goes from σ=0.15 at the 2 nd Order to σ=0.03 at the 3 rd Order. In addition, we can examine the variation of σ as a function of the 2 nd Order cutoff radius. This one is equal to 0.166, 0.164, 0.152, 0.151 and 0.150 when the cutoff radius varies from 8 a.u. to 16 a.u.. The introduction of the {100} shell reduces the standard deviation and leads to an improvement of the fit procedure.

In conclusion, an accurate description of the NTE in silicon can be achieved using a-TDEP, due to a good reproduction of its different parts (Bulk Modulus, thermodynamic Grüneisen parameter and specific heat). However, silicon does not show any strong variations of the phonon spectrum as a function of the temperature (in the [0 ; 400K] range). There is no (phase) stabilization or destabilization up to the melting. That is the reason why we show below two examples including explicit temperature effects that affect the phonon spectrum (and also the elasticity and the thermodynamic properties).

Thermal stabilization of the Zr bcc phase

In their original study, Hellman and coworkers [START_REF] Hellman | Lattice dynamics of anharmonic solids from first principles[END_REF] applied the TDEP method to the bcc phase of Zr. At 0 GPa, this structure is not stable at room temperature and all the calculations [START_REF] Wang | Firstprinciples calculations of phase transition, elastic modulus, and superconductivity under pressure for zirconium[END_REF][START_REF] Qian | Temperature effect on the phonon dispersion stability of zirconium by machine learning driven atomistic simulations[END_REF] performed at 0 K show few soft modes in the phonon spectrum around the ω point (corresponding to q = (2/3, 2/3, 2/3) along the H-P high-symmetry line) and the Γ point.

In the present work, AIMD simulations are performed with 128 atoms in the supercell, using the NVT ensemble, for two temperatures (300 and 1000 K) and two pressures (around 0 and 30 GPa) during several ps. Our calculations are performed by using the ABINIT code [START_REF] Gonze | Abinit: First-principles approach to material and nanosystem properties[END_REF] over hundreds to thousands of processors [START_REF] Bottin | Large-scale ab initio calculations based on three levels of parallelization[END_REF]. A projector augmented wave [START_REF] Blöchl | Projector augmented-wave method[END_REF][START_REF] Torrent | Implementation of the projector augmented-wave method in the abinit code: Application to the study of iron under pressure[END_REF] (PAW) atomic data is generated using the ATOMPAW code [START_REF] Holzwarth | A projector augmented wave (paw) code for electronic structure calculations, part i: atompaw for generating atom-centered functions[END_REF] In Fig. 8, we show the phonon spectra obtained by a-TDEP using only a 2 nd order expansion of the energy, without any higher order terms. The soft modes disclosed by a-TDEP at 300 K (responsible for the instability) are suppressed when the temperature increases. In addition, we also notice that the phonon spectrum obtained at 1000 K is in very good agreement with experiments [START_REF] Heiming | Phonon dispersion of the bcc phase of groupiv metals. ii. bcc zirconium, a model case of dynamical precursors of martensitic transitions[END_REF]. At last, we show in Fig. 8 the effect of pressure (30 GPa) on the phonon spectrum. In this case, the soft modes disappear which unveils the stabilization of the bcc phase at high pressure and room temperature, in line with experiments.

Failure of the QHA in α-U

Uranium is well-know to be the only element in the periodic table in which a charge density wave (CDW) state is observed at room pressure [START_REF] Lander | The solid-state properties of uranium a historical perspective and review[END_REF]. The phase transition at 50 K is accompanied by the softening of the longitudinal optic-like phonon mode in the [START_REF] Blöchl | Projector augmented-wave method[END_REF] direction [START_REF] Crummett | Lattice dynamics of α-uranium[END_REF][START_REF] Smith | Observation of a charge-density wave in α-u at low temperature[END_REF][START_REF] Bouchet | Lattice dynamics of α uranium[END_REF] and involves a doubling of the unit cell.

In the present work, AIMD simulations are performed with 96 atoms in a 4 × 2 × 3 α-U supercell, using the NVT ensemble, for three temperatures ( 300 and 900 K) during several ps. Our calculations are performed by using the ABINIT code [START_REF] Gonze | Abinit: First-principles approach to material and nanosystem properties[END_REF] over hundreds to thousands of processors [START_REF] Bottin | Large-scale ab initio calculations based on three levels of parallelization[END_REF]. A PAW [START_REF] Blöchl | Projector augmented-wave method[END_REF][START_REF] Torrent | Implementation of the projector augmented-wave method in the abinit code: Application to the study of iron under pressure[END_REF] atomic data is generated using the ATOMPAW code [START_REF] Holzwarth | A projector augmented wave (paw) code for electronic structure calculations, part i: atompaw for generating atom-centered functions[END_REF] in the XML format [103], considering 14 electrons in the valence and a cutoff radius of 2.85 bohr leading to a cutoff energy equal to 435 eV. Exchange and correlation energy was treated using the GGA-PBE functional [104]. In Fig. 9, we compare the temperature dependence of the longitudinal optic-like phonon mode obtained using the DFPT (at 0 K), with a-TDEP at 50, 300 and 900 K (using only a 2 nd order expansion of the energy, without any higher order terms) with experimental data at 50 and 300 K. We use the experimental lattice parameter at each temperature for both methods (for more computational details, please see Ref. [START_REF] Bouchet | Thermal evolution of vibrational properties of α-u[END_REF]).

With a-TDEP we capture the softening of the phonon mode. At 900 K the α-U structure is stable and as the temperature decreases from 900 to 50 K the phonon branch strongly softens with a well pronounced dip located near q = (1/2, 0, 0) along the [ζ00] high-symmetry line. The frequencies are in excellent agreement with neutron-scattering data at 300 and 30 K [START_REF] Crummett | Lattice dynamics of α-uranium[END_REF][START_REF] Smith | Observation of a charge-density wave in α-u at low temperature[END_REF].

On the other hand, the results obtained with DFPT in QHA are in contradiction with the experimental observations. The α-U structure is predicted to be unstable at volumes corresponding to 900 and 300 K, due to imaginary frequencies. By contrast, it becomes stable at volume corresponding to 50 K when the phase transition towards the CDW state takes place. This example clearly highlight the failure of the QHA when explicit temperature effects are present.

Conclusion

We have developed a post process code able to handle the anharmonic effects. This one is based on the TDEP method and is distributed in the ABINIT code under the GNU-GPL license. Many features can be computed using a-TDEP : phonon spectrum, elastic constants and moduli, Grüneisen parameters, thermal expansion, free energies... These capabilities have been checked on many materials (Si, U [START_REF] Bouchet | Thermal evolution of vibrational properties of α-u[END_REF][START_REF] Bouchet | High-temperature and high-pressure phase transitions in uranium[END_REF], Zr, Pu [START_REF] Dorado | Phonon spectra of plutonium at high temperatures[END_REF], MgO [START_REF] Bouchet | Ab initio calculations of the B1-B2 phase transition in MgO[END_REF], Ti [START_REF] Dewaele | High pressure-temperature phase diagram and equation of state of titanium[END_REF], Fe [START_REF] Dewaele | Iron under conditions close to the α -γ -triple point[END_REF]), in various thermodynamic situations (low/high temperature/pressure), and excellent agreements with experimental data are obtained everytime.

This paper paves the way to have a complete reference support for the user of a-TDEP. It includes the theory needed to understand the TDEP method and to compute thermodynamic quantities, the algorithms used in this a-TDEP implementation, numerous references to works applying this method but also to other similar approaches and finally, several examples showing the capabilities of this code.

The choice to implement the a-TDEP code in ABINIT is definitively valuable since this post-process actually benefits of all the developments from the ABINIT ecosystem : the symmetries coming from the lattice or magnetic space group, the treatment of the dipole-dipole interaction [START_REF] Gonze | Dynamical matrices, born effectives charges, dielectric permittivity tensors, and interatomic force constants from densityfunctional perturbation theory[END_REF], the input-output file formats (NetCDF, Yaml...), the definition of Monkhorst-Pack grids of q-points used to computed the vDOS or other thermodynamic quantities, the stochastic generation of thermalized configurations [START_REF] Kim | Nuclear quantum effect with pure anharmonicity and the anomalous thermal expansion of silicon[END_REF]... In the future, other interplay are expected, in particular the use of a-TDEP during the execution of ABINIT or MultiBinit [START_REF] Wojde L | First-principles model potentials for lattice-dynamical studies: general methodology and example of application to ferroic perovskite oxides[END_REF].

The present version of the a-TDEP code does not include several features. In particular, the 4 th order is not implemented [START_REF] Feng | Four-phonon scattering significantly reduces intrinsic thermal conductivity of solids[END_REF][START_REF] Yue | Metric for strong intrinsic fourth-order phonon anharmonicity[END_REF][START_REF] Ravichandran | Unified first-principles theory of thermal properties of insulators[END_REF], so the phonon lifetime, the frequency shift, the thermal conductivity, the spectral function and the anharmonic corrections to thermodynamic (free energy, heat capacity...) and elastic (sound speed, Grüneisen parameter...) quantities will be available in the next version.

Figure 1 :

 1 Figure 1: Sketch of interaction shells in 2-dimension. Left panel (second order) : six shells of pair interaction. The atoms belonging to each shell are connected through the crystal symmetries. Right panel (third order) : two kinds of triplet interactions.

  1 ...ip = Λp λp=1 Υ κ λp θ pλp with κ=1,...3 p with θ pλp the Λ p nonzero independent coefficients of the IFC ref p . The Υ κ λp matrix is built by gathering the column vectors orthogonal to the line vectors of Υ κ λp .

Figure 2 :

 2 Figure 2: qTep GUI for a-TDEP. First tab is for the supercell and MD trajectory, second one for the reference structure (unitcell) and the last one is for some available options.

Figure 3 :

 3 Figure 3: (Color online) Thermal expansion of silicon between 0 and 400 K. Experiments are from Carr et al. [73], Ibach et al. [74], Slack et al. [75], Lyon et al. [76] and Middelmann et al.[77].

29

 29

Figure 4 :

 4 Figure 4: (Color online) Variation of the elastic constants of silicon C 11 (black), C 12 (red) and the Bulk Modulus K T (blue) as a function of the 2 nd Order cutoff radius R cut2 . For each, a dotted line with the same colour mentions the experimental value.

Figure 5 :

 5 Figure 5: (Color online) Dispersion curves of the mode Grüneisen parameter of silicon computed using a-TDEP, for two cutoff radii (R cut3 =8 and 12 a.u.). The experimental data from Weinstein [80] are shown for comparison.

Figure 6 :

 6 Figure6: (Color online) Phonon spectra of silicon computed using a-TDEP for various values of 2 nd and 3 rd Order cutoff radii. The experimental data "Exp 1" are from Nilsson and Nelin[START_REF] Nilsson | Study of the homology between silicon and germanium by thermal-neutron spectrometry[END_REF] and "Exp 2" are from Kulda et al.[START_REF] Kulda | Inelastic-neutronscattering study of phonon eigenvectors and frequencies in si[END_REF].

Figure 7 :

 7 Figure 7: (Color online) Measure of the fit procedure : Forces of the model (TDEP) as a function of the AIMD ones for silicon. In the inset at the top left, we show the 2 nd Order, 3 rd Order and AIMD energies (they are all superimposed). In the inset at the bottom right, we show the energy differences (E 2nd T DEP -E AIM D ) and (E 3rd T DEP -E AIM D ).

Figure 8 :

 8 Figure8: (Color online) Phonon spectra of bcc zirconium computed using a-TDEP at two temperatures (300 and 1000 K) and around two pressures (0 and 30 GPa). Experiment data come from Ref.[START_REF] Heiming | Phonon dispersion of the bcc phase of groupiv metals. ii. bcc zirconium, a model case of dynamical precursors of martensitic transitions[END_REF].

Figure 9 :

 9 Figure 9: (Color online) Phonon frequencies of the longitudinal optic-like mode of α-U at 900, 300 and 50 K. Black lines : a-TDEP, red lines : DFPT. Experimental neutronscattering data [106] at room temperature are denoted by blue circles and by the dashed line at 30 K [107].

  θ 22 θ 23 θ 24 -θ 23 θ 22 θ 26 θ 25 θ 26 θ 22

  in the XML format [103], considering 12 electrons in the valence and a cutoff radius of

		6
	Frequency (THz)	2 4	300 K -0 GPa -30 GPa
		0	1000 K -0 GPa -30 GPa
			Exp
		Γ -2	H	ω P	Γ	N

or 16 a.u. second one. In all cases, the agreement with experiments[START_REF] Weinstein | Raman scattering and phonon dispersion in si and gap at very high pressure[END_REF] is very good. In the same time, this affects the thermodynamic Grüneisen parameter which goes from 0.34 (R cut3 =8 a.u.) to 0.42 (R cut3 =12 a.u.) and 0.45 (R cut3 =16 a.u.).The third contribution to the thermal expansion is the specific heat C V (V, T ), built using the phonon spectrum ω(V, T ). As displayed in Fig.6, the variation of the phonon dispersion curves as a function of the 2 nd order

2.2 bohr leading to a cutoff energy equal to 400 eV. Exchange and correlation energy was treated with the generalized gradient approximation (GGA) using the Perdew-Burke-Ernzerhof (PBE) functional[104]. As for Si, a careful treatment of the electronic density integration has been carried out with the use of a (2×2×2) k-point mesh.

Acknowledgments

We would like to thank Olle Hellman, Matthieu Verstraete and Xavier Gonze for all the fruitful discussions we had, and Antoine Levitt for assisting us in the development and implementation of computational methods. We are grateful to the two unknown referees for their remarks which significantly improve the quality of the manuscript.


Choron Guillaume

Lesage François-Xavier

Picy Laurence

Courtet Philippe

Olie Emilie

Emilie Olie
email: e-olie@chu-montpellier.fr

Evaluation somatique de cent patients hospitalisés en unité de crise psychiatrique : une étude observationnelle rétrospective Somatic assessment of one hundred inpatients in a psychiatric crisis unit: a retrospective observational study

Keywords: Mental disorders, Mortality, Comorbidity, Health status, Cardiovascular risk 43 Troubles mentaux, Mortalité, Comorbidité, Etat de santé, Risque cardiovasculaire

Objectives: Compared to the general population, psychiatric patients are more exposed to physical illness but have reduced access to care. Methods: We conducted a descriptive study in Montpellier between November 2011,2nd and December 2012, 21st. Every Wednesday and Friday, the last two inpatients admitted in the psychiatric post-emergency unit of the University-Hospital of Montpellier were assessed by a general practitioner and included in the study. This unit takes care of suicidal patients suffering from mood and / or personality disorders. The general practitioner assessed lifetime somatic history, current somatic comorbidities and medical follow-up for non-psychiatric issues. Results: One hundred patients were included. The sample consisted of 81% women with a mean age of 43 years (18-79 years). The majority of patients had a lifetime history of somatic disease (96%) and was followed by a general practitioner (99%). Six patients had no met general practitioner for at least one year. Dyslipidemia was reported in 32 patients, among those only one patient was on lipid-lowering drug (96.88%). Among patients with impaired dental condition, 29 (55.77%) had not met a dentist for at least one year. Among those with impaired near visual acuity and/or impaired distance visual acuity, 19 (65.52%) had not met an ophthalmologist for at least one year. Conclusion: Although detected, somatic comorbidities seem outsourced in psychiatric patients. Greater awareness among different health professionals about the medical care of such patients could improve healthcare management and life expectancy.

RESUME

Intérêt : En comparaison de la population générale, les personnes souffrant de troubles psychiatriques sont en moins bonne santé et ont un accès aux soins réduits. Méthodes : Nous avons réalisé une étude descriptive au CHU de Montpellier entre le 2 novembre 2011 et le 21 décembre 2012. Chaque mercredi et vendredi, les deux derniers patients hospitalisés dans le service de post-urgences psychiatriques du CHRU de Montpellier étaient évalués par un médecin généraliste et inclus dans l'étude. Ce service accueille principalement des patients en crise suicidaire dans le cadre d'un trouble de l'humeur et/ou d'un trouble de personnalité. L'évaluation informait sur les antécédents somatiques, les comorbidités somatiques actuelles et le suivi médical pour des raisons non psychiatriques. Résultats : Cents patients ont été inclus. La population était principalement composée de femmes (81%) et l'âge moyen était de

INTRODUCTION

Compared to the general population, psychiatric patients are more exposed to physical illness because of: (i) general risk factors (genetic, environmental, occupational…), (ii) psychotropic side-effects, (iii) risk factors related to the patient psychopathology and (iv) risk factors related to psychiatric comorbidities (addiction…). Hence, psychiatric patients have a higher premature mortality rate, an unhealthy lifestyle and poorer access to health care services [START_REF] Harris | Excess mortality of mental disorder[END_REF][START_REF] Colton | Congruencies in Increased Mortality Rates, Years of Potential Life Lost, and Causes of Death Among Public Mental Health Clients in Eight States[END_REF][START_REF]Tilgangen til somatiske helsetjenester blant psykisk syke[END_REF]. The main risk identified is cardiovascular risk which increases the overall mortality of psychiatric patients [START_REF] Charrel | Mortality of people suffering from mental illness: a study of a cohort of patients hospitalised in psychiatry in the North of France[END_REF]. We need to have a comprehensive approach and enhance the dialogue between general practitioner and other specialists. General practitioners are the first doctors consulted by patient suffering from psychological distress: more than 60% patient suffer from psychological distress would see a general practitioner versus 10% who would see a psychiatrist [START_REF] Lépine | Depression in the community: the first pan-European study DEPRES (Depression Research in European Society)[END_REF]. Furthermore, general practitioners are often the gateway to the health care system for psychiatric patients [START_REF] Dezetter | Use of general practitioners versus mental health professionals in six European countries: the decisive role of the organization of mental health-care systems[END_REF]. Nevertheless, studies showed that the quality of communication between psychiatrists and general practitioners needs improvement [START_REF] Colton | Congruencies in Increased Mortality Rates, Years of Potential Life Lost, and Causes of Death Among Public Mental Health Clients in Eight States[END_REF][START_REF] Cohidon | La non-communication entre la médecine générale et le secteur de santé mentale[END_REF][START_REF] Aubin-Auger | La consultation du patient à risque suicidaire en médecine générale Généralistes et psychiatres : une relation compliquée[END_REF]. Fleury et al. showed that collaboration between general practitioner and psychiatrists was a key factor for accurate management of mental disorders whereas limited access to specialist was impediment [START_REF] Fleury | General practitioners' management of mental disorders: A rewarding practice with considerable obstacles[END_REF]. The medical care of patient suffering from severe mental disorders has been the subject in France of « Haute Autorité de Santé » (HAS) recommendations, recently proposed by the FFP-CNPP (Fédération Française de Psychiatrie -Conseil National Professionnel de Psychiatrie [START_REF] Elm | The Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) statement: guidelines for reporting observational studies[END_REF] . Thus, we aimed at describing the general health status of inpatients in a psychiatric post-emergency unit in order to identify the elements of monitoring that could be improved.

METHODS

A random sample of 100 inpatients were analyzed retrospectively between November 2011, 2nd and December 2012, 21st. The department of Psychiatric Emergency and Post-acute Care mainly takes care of patients in suicidal crisis in the context of a mood disorder and / or borderline personality disorder.

Two days a week (Wednesday and Friday), a standardized assessment of the two last admitted inpatients was performed by a general practitioner. This assessment included a full clinical examination that provided information on (i) somatic personal history, (ii) current somatic comorbidities (including cardiovascular comorbidities), (iii) follow-up by a general practitioner (frequency, date of last visit) and (iv) follow-up by a specialist (gynecologist, ophthalmologist and dentist). A biological sample was performed in order to detect metabolic disorders.

Qualitative variables were described according to frequency and percentage. Quantitative variables were described by the mean and standard deviation or median and interquartile range when graphical distribution was normal or not respectively. Descriptive analysis of the data was be performed using the R software for Windows version 3.3.0. STROBE guidelines were followed (11).

RESULTS

In total, one hundred inpatients were included (Table1). Most of the patients were women (81%), middle-aged (mean age = 41 years; SD = 14.75), single (60%) and had children (65%). Half of the patients were working full-time job (50%) and 17% were student. The mean age at onset of psychiatric disorder was 18.5 years (SD = 9.5) while the mean age at first medical contact for psychiatric reasons was 24.5 years (SD = 10.3). Only half of the patients had a consultation with a psychiatrist within the year following psychiatric disorder onset. The mean-age at first psychiatric hospitalization was 30 years (SD = 12.66).

The most frequent comorbidities were cardiovascular, including hypertension (17%), orthostatic hypotension (19%), cardiac murmur (11%), dyspnea (16%) and lower extremity oedema (9%). Half of patients (45%) were smokers. Two third of patients did not practice physical activity. Thirteen patients were underweight (body mass index < 18.5 kg/m²), 20 were overweight (body mass index between 25.0 and 29.9 kg/m²) and 9 were obese (body mass index > 30 kg/m²). One third of patients (36%) had a high waist-to-hip ratio (> 0.85 for females or > 1 for males). Over one third of patients (39%) had a high waist size (> 80 cm for females or > 94 cm for males). A sixth of patients (14%) had an excessive consumption of alcohol defined by > 14 units per week for females and > 21 units per week for males. Among eight patients with low blood sugar, 3 patients were on anti-diabetic drugs. One fifth of patients had hyperglycemia and/or were treated with antidiabetic drugs. One third of patients (32%) had dyslipidemia (hypercholesterolemia or hypertriglyceridemia). Nine patients were on lipid-lowering drugs among those one had uncontrolled hyperlipidemia. Nineteen patients had a metabolic syndrome according to the International Diabetes Federation criteria [START_REF] Alberti | The metabolic syndrome-a new worldwide definition[END_REF][START_REF] Boursier | Le syndrome métabolique[END_REF] .

Most of patients (96%) had a health insurance. Almost all patients (99%) reported having a regular general practitioner. But six patients had not visited their general practitioner for more than one year.

Oral examination revealed moderate to severe alterations in 52 patients, among those, half had not visited dentist within the past year. Five patients have never been to a dentist. Two thirds of patient suffering from low visual acuity (65.52%) had not visited an ophthalmologist within the past year.

Among the 81 females, 14 had not visited a gynecologist for more than 3 years. Five had never been to a gynecologist. Half of non-menopausal women did not use contraception (57.81%).

DISCUSSION

We show (i) almost all patients had a regular general practitioner, (ii) prescription rate of lipid-lowering drug was low despite high rate of dyslipidemia, (iii) visits to specialists was less frequent in psychiatric sample than in the general population.

In our sample, the prevalence of metabolic syndrome was higher than in the general population [START_REF] Junquero | Syndrome métabolique : quelle définition pour quel(s) traitement(s) ? ms[END_REF]. Moreover, we may have underestimated the frequency of metabolic syndrome as thirteen patients did not have a blood test. Presence of metabolic syndrome is associated with an increase mortality in psychiatric patients [START_REF] Łopuszańska | Mental illness and metabolic syndrome -a literature review[END_REF]. Intriguingly, most of patients with dyslipidemia did not receive adequate pharmacological treatments. This result is consistent with a previous observation done in Fondamental Advanced Centers of Expertise in Bipolar Disorders (FACE-BD) (16). It suggests that metabolic disorders are enough screened in primary care even if psychiatric disorder are at high metabolic risk. Recommendations [START_REF] Saravane | Élaboration de recommandations pour le suivi somatique des patients atteints de pathologie mentale sévère[END_REF] highlight the necessity to realize regular clinical, biological and electrocardiographical examinations in subjects suffering from psychiatric disorder. Furthermore, many psychotropic drugs participate to metabolic disorders, imposing annual biological monitoring. Several studies showed general practitioners have difficulties to take care of patients suffering from psychiatric disorders. Emergency care practitioners feel insufficiently trained to manage suicidal patients [START_REF] Betz | Knowledge, Attitudes and Practices of Emergency Department Providers in the Care of Suicidal Patients[END_REF]. Some general practitioners declare to feel so much uncomfortable with the psychiatric patients that they "avoid" them. Some others general practitioners would doubt about the quality of delivered medical care. These difficulties, which may have prognostic consequences, were well reported in a report by French Senat in 2009: "Delayed diagnosis of psychiatric disorders and overconsumption of psychotropic drugs are partly due to the lack of psychiatry information delivered to the general practitioners" [START_REF]La psychiatrie en France : de la stigmatisation à la médecine de pointe[END_REF]. Somatic management of psychiatric patients needs to be integrated in training of general practitioner. On the contrary, some psychiatrists underestimate the cardiovascular risks of their patients. European Society of Cardiology listed the recommendations for the lipid-lowering drug treatments of patients suffering from mental disorders: (i) major psychiatric disorders are modifiers for estimating total cardiovascular risk, (ii) the management of total cardiovascular risk in patients with a psychiatric disorder is not different from what is recommended in patient at high/very high cardiovascular risk, (iii) in patients with psychiatric disorders particular attention has to be paid to adherence to lifestyle changes and compliance with drug treatment [START_REF] Reiner | ESC/EAS Guidelines for the management of dyslipidaemias: The Task Force for the management of dyslipidaemias of the European Society of Cardiology (ESC) and the European Atherosclerosis Society (EAS)[END_REF]. The European Society of Cardiology proposed different interventions according to cardiovascular status and level of low-density lipoprotein cholesterol [START_REF] Reiner | ESC/EAS Guidelines for the management of dyslipidaemias: The Task Force for the management of dyslipidaemias of the European Society of Cardiology (ESC) and the European Atherosclerosis Society (EAS)[END_REF]. Systemic prescription of statins should be considered in psychiatric patients with dyslipidemia to manage cardiovascular risk [START_REF]The Effect of Concomitant Treatment With SSRIs and Statins: A Population-Based Study[END_REF].

The rate of visits to general practitioner is similar in our sample and in general population [START_REF]The Effect of Concomitant Treatment With SSRIs and Statins: A Population-Based Study[END_REF], but psychiatric patients had a less frequent access to specialists (gynecologists, ophthalmologist and dentists) than general population. For example, 29% patients had not visited a dentist within the two last years vs 18% in the general population. We propose two explanations: depressive symptomatology of depression and precariousness [START_REF]Le score Épices : un score individuel de précarité[END_REF].

Most of the patients were women (81%), middle-aged from the South of France. This selection bias impact the generalizability of our results.

CONCLUSION

Although detected, somatic comorbidities seem outsourced in psychiatric patients. Greater awareness among different health professionals about the medical care of psychiatric patients could improve healthcare management and life expectancy.

Conflit d'intérêt : aucun