Keywords: porous carbon, carbon nanostructures, simulation, adsorption theory

  L'archive ouverte pluridisciplinaire

Editorial

The objective of this Special Issue is to present various modern aspects of computer simulation applied to the study of porous and nanotextured carbon materials including zeolites, nanotubes, soot particles but also more conventional disordered (activated or not) carbons. Carbon materials are at the heart of many high-tech industrial processes in electronics, aerospace technology, energy storage and filtration but are also aerosol air pollutants when emitted by cars or airplanes. It is therefore essential to acquire knowledge of the texture of these various carbon forms and related properties at the nanometric scale or even below. Carbons materials are known as very efficient adsorbents due to their high specific surface area and porosity. However the concept of pore in the real carbon materials is far from being obvious. Even in the case of carbon nanotubes that have individually a well-defined cylindrical geometry, the way by which the porosity distributes among a bundle of tubes is not yet clear. This is critical in electronic or field emission applications where dense carpets of well-aligned tubes are required. The situation gets more difficult when considering common porous carbons obtained from the thermal decomposition of organic precursors; saccharose cokes for instance. These materials all have a degree of disorder at various length scales and exhibit local nanostructures.

For characterization purposes, experimentalists use a combination of techniques that often includes low-temperature gas adsorption and transmission electron microscopy (TEM). Information from adsorption-based methods is gained assuming that the pore network can be reduced to a distribution of pores (PSD) of a given characteristic size and geometry (slit, cylinder). By contrast, TEM is a direct measure of some aspects of the material texture but 3D information such as the pore morphology (shape) and topology (the way pores distribute in space) is hard to obtain since TEM images are projections of 3D structures under an electron beam. Since the early 70's, carboneous materials TEM data have been analysed in terms of Basic Structural Units (graphene plane segments) of various sizes. These structural features seem common to a lot of carbon phases. This is the reason why theoreticians have been so keen in modelling porous carbons as an assembly of slit pores of various sizes with graphitic walls of various thicknesses. Starting in the 80's, this model was a great help in elucidating fundamental properties of confined fluids and related phenomena as a function of the pore size: micropore volume filling mechanism to capillary condensation in nanopores. These are still the subject of intense research. An objective of this Special Issue was to bring together experts in different approaches, with a view to understanding the relation between the different methods in studying and predicting adsorption properties. Several questions can be addressed. What are the strengths and weaknesses of the various available approaches and models? What are the applications for which a given approach is most suited? The reader will find in the present Special Issue several simulation and theoretical fundamental studies describing the adsorption properties of polar, non polar fluids and mixtures in slit pores or graphite surfaces [Do. et al ; Quirke et al ; de Weireld et al].

A few years ago, the first simulation studies of carbon porous networks not based on the slit pore model but aiming at giving some account of the real pore morphology and topology were published showing that atomistic level simulation methods can be very helpful in The last paper of this Special Issue tackles the question of the validity of one of most commonly used adsorption-based characterization method for pore PSD determination (namely the Density Functional Theory) that is available on most adsorption commercial machines [Roussel et al]

  Atomistic Simulation and Theory of Nanoporous Carbons and Nanostructures

  of disordered materials. However, in many cases, these also are beset with difficulties: the improvement of such atomistic approaches over the slit pore model in predicting adsorption and separation properties is not yet clear. The time has come to tackle this issue. The reader will find in the present Special Issue several papers on atomistic realistic models of porous carbon[Do et al ; Bhatia et al ; Biggs et al] along with some simulation studies of adsorption and freezing of fluids in various realistic models of porous carbons ranging from activated saccharose cokes[Coasne et al] to soot particles[Picaud et al].

Page 1 of 2 http://mc.manuscriptcentral.com/tandf/jenmol

	1
	2
	3
	4
	5
	6
	7
	8
	9
	10
	11
	12
	13
	14
	15
	16
	17
	18
	19
	20
	21
	22
	23
	24
	25
	26
	27
	28
	29
	30
	31
	32
	33
	34
	35
	36
	37
	38
	39
	40
	41
	42
	43
	44
	45
	46
	47
	48
	49
	50
	51
	52
	53
	54
	55
	56
	57
	58
	59

Page 2 of 2 http://mc.manuscriptcentral.com/tandf/jenmol

  and routinely used in industry and research laboratories.

	1
	2
	3
	4
	5
	6
	7
	8
	9
	10
	11
	12
	13
	14
	15
	16
	17
	18
	19
	20	Roland J.-M. Pellenq
	21
	22
	23
	24
	25
	26
	27
	28
	29
	30
	31
	32
	33
	34
	35
	36
	37
	38
	39
	40
	41
	42
	43
	44
	45
	46
	47
	48
	49
	50
	51
	52
	53
	54
	55
	56
	57
	58
	59
	60


Frederic Affouard

Patrice Bordat

Adrien Lerbret

Marc Descamps

Slow dynamics in glass-forming materials

Keywords: orientational disordered crystals, glassforming liquids, fragility, liquid-glass transition, Mode-Coupling Theory orientational disordered crystals, glassforming liquids, fragility, liquid-glass transition, Mode-Coupling Theory, Coupling Model, anharmonicity. Molecular Simulation

published or not. The documents may come

Introduction

It is well known that very different materials such as: silica, low-molecular-weight liquids (carbohydrates, alcohols), polymers or even proteins are able to exhibit a very intriguing feature called glass transition (1). This latter is characterized by an extraordinary decrease of several orders of magnitude of the mobility in a narrow temperature range without significant structural changes.

Understanding of the mechanisms leading to the glass transition has been a subject of intense research in recent years and it is still a matter of numerous theoretical and experimental studies (1).

Much attention has been recently devoted to the high temperature pico-nanosecond regime (ps-ns) of liquid glass-formers. A fundamental question concerns the microscopic description of the cooperative mechanisms which develop over a temperature range in which dynamics start being both relatively well described by the predictions of the mode coupling theory (MCT) and influenced by the potential energy landscape (2; 3). So far, MCT (4) is the only theory which provides a microscopic description of supercooled atomic liquids. The intrinsic basis of MCT states that the behavior of any time-dependent correlators describing the dynamics of the system is only controlled by its static density correlator S(Q) and predicts the existence of a critical temperature T c corresponding to an ergodic to non-ergodic transition. Scaling properties predicted by MCT approaching T c have been successfully validated in numerous experiments (5; 6; 7; 8) and molecular dynamics (MD) simulations [START_REF] Schilling | Theories of the structural glass transition[END_REF].

It is now well accepted that T c also marks the crossover to the "landscape dominated" regime, as advocated long time ago by Goldstein [START_REF] Goldstein | [END_REF], where dynamics can be described as thermally activated hoppings between inherent structures (minima of the potential energy). It has been demonstrated that an other remarkable crossover temperature T A can be identified above T c which coincides with the onset of non-exponential relaxation i.e the "landscape influenced" regime (11; 3).

The aim of this paper is to present and discuss results obtained from MD simulations on slow dynamics The idealized version (4) of this theory predicts a two-step relaxation scenario (fast β, slow α) of all the time dependent correlators φ(q, t). MCT particularly states the following points: At short time, φ(q, t) decays to a plateau value, the so-called nonergodicity parameter classically noted f c q . The dynamical regime associated with this plateau is called β. It is centered around the rescaling time t σ which is given by t σ = t 0 |σ| -1/2a where t 0 is a characteristic microscopic time and σ is proportional to |T -T c |. Above T c , the late β regime or the early α relaxation is described by the following power law (going beyond first order): φ(q, t) = f c q -h (1) q .(t/τ ) b + h (2) q .(t/τ ) 2b

where the first two terms correspond to the classical von Schweidler law, the last term is a second order correction and τ = t 0 |σ| -γ with γ = 1/2a + 1/2b. At long time, MCT also predicts that the previous fast regime is followed by a slow relaxation, called α, with the characteristic time τ . Furthermore, parameters a and b are temperature and q independent and related via Γ

2 (1 -a)/Γ(1 -2a) = Γ 2 (1 + b)/Γ(1 + 2b)
where Γ(x) is the gamma function.

Glassy Crystals

Most of the glassformers are molecular liquids and accordingly possess both translational and rotational degrees of freedom (TDOF and ODOF). Approaching the glass transition, the coupling between TDOF and ODOF is far from being understood and the precise role of the ODOF remains unclear as demonstrated by the so-called low-temperature translation-rotation paradox (12; 13). It is now well established that the phenomenology of glass formation is also displayed by some partially disordered molecular crystals i.e. plastic crystals in which the average position of the centers of mass are ordered on a lattice while the orientations are dynamically disordered (14; 15). Indeed, some of them, called glassy crystals [START_REF] Suga | [END_REF], such as cyanoadamantane (17), ethanol (18; 19), cyclooctanol (20) or difluorotetrachloroethane (21) can be deeply supercooled and present many properties characteristic of the conventional molecular liquid glasses such as a step in the specific heat at the glass transition or a non-Arrhenius behavior of the relaxation times. Those compounds has led recently to a strong interest since they offer valuable possibilities to focus mainly on the role of the orientational degrees of freedom (ODOF) whose the importance during the glass formation remains a matter of debate at present. It has to be noticed that glassy crystals in opposition to the so called "orientational glasses" like mixed crystals Ar 1-x (N 2 ) x or (KCN) x (KBr) 1-x are not frustrated by quenched disorder externally imposed by a dilution (22). MCT has been initially developed to provide a microscopic understanding of simple atomic liquids. However, general predictions have also been successfully validated for molecular compounds. Recently, some extensions of MCT called molecular mode coupling theory (MMCT) [START_REF] Schilling | Theories of the structural glass transition[END_REF] have been proposed to take ODOF into account: one diatomic probe molecule in an atomic liquid (23), liquids made of linear molecules (24) or water (25) to cite only a few. The authors have particularly shown that some of the basic predictions of MCT still hold owing TDOF/ODOF coupling. As revealed by recent MD investigations (26; 27) performed on orthoterphenyl (OTP), coupling of the rotational dynamics to the center-of-mass motion can be complex. No microscopic theory has been developed for plastic crystals so far. However, in (28; 29; 30) In the following, we investigated three models of plastic crystals possessing different kind of crystalline and molecular symmetries and different type of reorientational dynamics. The present numerical simulation investigations aim to demonstrate that plastic crystals whose the dynamics are almost completely controlled by rotation shares some common dynamical features with liquid glass-formers in the ps-ns regime where some precursor features (non-exponentiality or non-Arrhenius behavior) of the glass transition emerge.

Models

3.1.1

Chloroadamantane. Chloroadamantane C 10 H 15 Cl (CLADM) is a rather huge molecule which belongs to the substituted adamantane family. It shows a plastic phase structure isomorphous to cyanoadamantane, but the chloroadamantane molecule possesses a smaller substitute and a faster dynamics well adapted for MD simulation investigations. Rotational motions are also suspected to change in nature as it is reported from earlier incoherent quasielastic neutron scattering experiments (31). Chloroadamantane undergoes at 244 K a first order transition from an ordered monoclinic structure to a rotator phase with face-centered-cubic (fcc) symmetry (32). The plastic-liquid transition occurs at 442 K (32).

The simulated system is composed of rigid linear molecules with two sites: one chlorine atom (noted Cl) and one super atom (noted Adm) that models the adamantane part C 10 H 15 . The distance of both pointlike atoms is 3.3 Å and the molecule's moment of inertia is 302.733 amu. Å2 . Rigid body equations with quaternions were used for the the orientational variables.Molecular dynamics calculations were performed on a system of N = 256 molecules (4 × 4 × 4 fcc crystalline cells) interacting through a Lennard-Jones short range site-site potential of the form

v(r) = 4 ((σ/r) p -(σ/r) q ) ( 2
)
where r is the distance between two different sites. The parameters , σ, p and q are specified in table 1.

A complete description of a similar model used previously to simulate cyanoadamantane using the Monte Carlo technique is given in (33; 34), so we give here only the essential details. The parameters of the Cl-Cl, C-C and H-H potentials were taken from potential energy calculations of chlorine atoms in molecular crystals (35). The parameters for Adm interactions have been determined from the atom-atom interactions between atomic sites of the adamantane part of the molecules (36). First, the averaged energies over all molecular orientations were calculated for different distances of centers of mass of the adamantane part of molecules. Then, they were fitted into the Lennard-Jones form (seei table 1). The parameters for interactions between dissimilar sites have been determined by the Berthelot combining rules. The 12-6 Lennard-Jones functions were applied for all dissimilar interactions.

The chloroadamantane molecule possesses a relatively large dipolar moment µ (2.39 Debyes) which is parallel to the molecular axis. The electrostatic interactions were handled by the Ewald method with two partial charges (q = ±0.151e) localized on both sites. This simple charge distribution has been constructed to make the total molecular charge equal to zero and to reproduce the total dipole moment. Newton's equations of motion were solved with a time step of ∆t = 5 fs. We worked in the NPT (constant number of molecules, temperature and pressure) statistical ensemble using the Melchionna modification of the Hoover algorithm with a thermostat and barostat relaxation time of 1 ps (37). Periodic boundaries conditions were used. MD simulations have been performed at 29 different temperatures from T = 220 to 500 K for a sample corresponding to the crystalline fcc rotator phase.

It should be noted that the CLADM model is particularly oversimplified. The main reason of this simplification is due to the large dipole moment of the CLADM molecule which requires the calculation of the electrostatic interactions unlikely to the Norbornylene or the difluorotetrachloroethane models (see below). However, it has been precisely constructed to keep the main experimental features and to speed up simulations in order to reach the nanosecond regime. Several comparisons of the results obtained from the CLADM model with experimental data have been published in (28). In this study, from calculations of the enthalpy, it was shown that this model captures the anomaly found in calorimetric experiments by Oguni et Al. (38) close to the Arrhenius to Non-Arrhenius deviation reported in the present manuscript. Moreover, both high-and low-temperature apparent activation energies of the relaxational time as defined in the present manuscript were found to be in good agreement with NMR experimental results. In (39), it is also reported that the shrinking of the lattice constant found in good agreement with X-ray experiments leads to an increase of the steric hindrance with decreasing temperature and in turn to a slowing down of the orientational dynamics.

3.1.2

Norbornylene. Norbornylene (NBE) C 7 H 10 displays one plastically crystalline over the temperature range 129 K to 320 K [START_REF] Westrum | Molecular Dynamics and Structure of Solids[END_REF]. The plastic crystal structure is hexagonal close-packed (hcp) [START_REF] Strange | [END_REF]. The structure of the brittle phase is reported to be monoclinic P 2 1 /c from ab-initio calculation (42). At high temperature, an isotropic reorientational motion is found (43) whose relaxational time follows an Arrhenius law. No discontinuity for this process is found at the plastic to liquid transition. Below 160 K, an abrupt dynamical change is mentioned with anisotropy setting.

The simulated sample consists of N = 686 molecules located on an average hcp lattice corresponding to a sample of 7 × 7 × 7 unit cells. Each molecules C 7 H 10 is described by its 17 atoms and considered as rigid.Rigid body equations with quaternions were used for the the orientational variables.Due to the presence of a carbon-carbon double bond, NBE possesses a small dipolar moment which was determined to be 0.32 D. As already made by Min et al. (42) in their ab initio determination of the NBE low temperature crystalline structure, the electrostatic contributions have been neglected in our simulation. Therefore, we have just considered a Buckingham atom-atom potential whose coefficients are given in table 2. Despite this approximation, both structure and dynamics were found in good agreement with experimental results. Newton's equations of motion were solved with a time step of ∆t = 5 fs. We worked in the NPT statistical ensemble with periodic boundary conditions at constant atmospheric pressure. We used Berendsen barostat and thermostat with relaxation times of 2.0 ps and 0.2 ps for controling the pressure and the temperature respectively (44). MD simulations have been carried out at 28 different temperatures from T = 125 to 320 K for a sample corresponding to the experimental crystalline hcp rotator phase. Below 125 K, the system falls out of equilibrium over the MD time window and becomes non ergodic. (45; 46; 47). DFTCE has been experimentally widely studied and presents a rich variety of unusual properties. DFTCE exhibits a glass transition of the overall rotation of the molecules at T g = 86 K (48). Two additional heat anomalies are found at 60 and 130 K and associated respectively to a sub-T g β process and the freezing of the transformation between trans and gauche conformation of the molecule. Changes concerning the nature of dynamics in this system has been reported from NMR experiments (49), Brillouin and dielectric spectroscopy (21). In this latter investigation, a crossover from individual to correlated rotational motions was particularly suggested.

MD calculations were performed on a system of N = 686 (7 × 7 × 7 bcc crystalline cells) molecules. Each DFTCE molecule is described by its 8 atoms and considered as a rigid unit. Rigid body equations with quaternions were used for the the orientational variables.They interact through a Buckingham short range atom-atom potential (see parameters in table 3) and the electrostatic contributions have been neglected since the DFTCE molecule possesses a weak dipolar moment. No structural or dynamical change has been found for a system where electrostatic interactions are taken into account. Newton's equations of motion were solved with a time step of ∆t = 5 fs. We worked in the NPT statistical ensemble with periodic boundary conditions where the simulation box is allowed to change in size and shape using the Melchionna modification of the Hoover algorithm (37) with a thermostat and barostat relaxation time of 1 ps. MD simulations were done for a sample corresponding to the DFTCE orientationally disordered phase at 14 different temperatures from T = 130 to 260 K in steps of 10 K. It should be mentioned that our very simple DFTCE model allows us to perform very long MD runs of about 50 ns. Single-molecule reorientational dynamics can be investigated from the self angular correlation functions C l=1,2 defined as:

C l (t) = 1 N i=1,N P l ( u i (t). u i (0)) (3)
where P l is the l-order Legendre polynomial and u i a unit vector parallel to the dipole moment of each molecule i in the case of CLADM and NBE or directed along the C 3 -symmetry carbon-carbon molecular axis of molecule i in the case of the DFTCE molecules. From C l=1,2 , we can obtain clarification on the nature of dynamical changes and anisotropy of the molecular motions. Experimentally, C 1 (t) can be directly measured in dielectric relaxation and C 2 (t) in Raman scattering [START_REF] Sherwood | The Plastically Crystalline State[END_REF]. C l=1,2 can also be related to the informations obtained from NMR relaxation measurements.

Figure 1 shows the C 2 (t) time correlation functions for all models and investigated temperatures. Clearly, at least for the lowest temperatures, a two step relaxation behavior as predicted by MCT and already observed in supercooled molecular liquids (45; 50) is shown in our simulations. At very short times, we observe a fast decay corresponding to a microscopic regime which is modulated by vibrational motions. This latter behavior is followed by a plateau-like region, which becomes more pronounced at the lowest temperatures. This feature is identified with the β regime. At the end, the last step is a decay to zero corresponding to the α process. C 1 (t) time correlation functions, not displayed in the present paper, show a very similar behavior to the C 2 (t) functions. When lowering the temperature, an intermediate plateau region emerges which proves the presence of an orientational caging between neighboring molecules (see Fig. 2). This is the rotational analogue of the translational cage effect observed in liquids (see Ref. (51)). This transient regime is followed by the α process which can be associated in plastic crystals with large tumbling motion between preferred molecular orientations as displayed in Fig. 2 for CLADM. A consequence is that the tumbling of one molecule is allowed only if an orientational rearrangement of its local neighbors occurs. This cooperative motion is clearly displayed using projections of the individual dipolar moments on one crystallographic plane.

The τ 1 and τ 2 relaxational times are defined as the time it takes for their respective time correlation functions to decay e -1 of their initial values. At high temperature, for all correlators, τ (T ) can be well fitted with an Arrhenius law τ (T ) = τ 0 exp(E 0 /T ) where the parameters (τ 0 , E 0 ) allow us to calculate the temperature dependent activation energy E(T ) = T.log(τ /τ 0 ).

Figure 3 displays E(T )/E 0 as function of the temperature. Clearly, at a remarkable temperature identified to T A , all activation energies start diverging simultaneously from unity corresponding to the Arrhenian high temperature behavior. This is an indication of the existence of a dynamical crossover from a free diffusion regime to the so-called landscape-influenced regime occurring at T A corresponding to a transition from quasi-free rotational diffusion of the molecular dipoles to activated geared tumblings as seen in Fig. 2.

The nature of the rotational motions involved in this dynamical change is given by the ratio τ 1 /τ 2 in Fig. 4.

Upon lowering the temperature, we observe that τ 1 /τ 2 clearly deviates from 3 and approaches 1. For uncorrelated small angular steps motion, it can be shown that the orientation correlation functions follow an exponential decay C l = exp[-l(l + 1)D r t] where D r is the rotational diffusion coefficient. τ 1 /τ 2 = 3 is therefore expected for free small-step rotation diffusion type of motion while τ 1 /τ 2 = 1 is classically associated to activated jump-like motion (52). These results prove the existence of a plastic-plastic transition which can be interpreted as the rotational analogue of the Goldstein crossing temperature seen in liquids.

Glycerol molecular liquid

We now focus on the dynamical properties of molecular liquids which present a slowing down upon cooling. Among them, we have studied glycerol for its key role in bioprotection. Indeed, it is well known that disaccharide/water binary solutions are able to protect biomolecules against dessication, dehydration and this ability seems to be strengthened by diluting a small amount of glycerol in these sugar/water solutions. However, the microscopic mechanism responsible for the biopreservation are still unclear and several hypotheses non fully satisfactory have been proposed (53; 54; 55; 56). We have realized a comparative study of the disaccharide/water solutions (57; 58) by looking at both the structure and the dynamics in order to extract characteristic parameters playing a key role in the bioprotection mechanisms. The next step is the study of ternary mixtures first implying the development of a good model of the glycerol liquid to reproduce the main dynamical features observed experimentally and presented below.

Computational details

The system is composed of 108 flexible glycerol molecules modeled using the all-atom AMBER force field (59) employed by Chelli et al. (60). This model produces a structure in good agreement with the experimental one but the dynamical properties can be improved. Then, the atomic charges and the Lennard-Jones parameters of the hydroxyl groups have been reparametrized (61) in order to reproduce more realistically the diffusion coefficient D (62) and the characteristic times τ Q obtained from Neutron Spin Echo experiments (63). Electrostatic interactions were handled by the reaction-field method with ε RF = 40 close to the values found in the literature [START_REF] David | CRC Handbook of Chemistry and Physics[END_REF]. The reaction field method has been compared to the Ewald summation on glycerol liquid and neither structural nor dynamical changes have been noticed. As glycerol molecules are entirely flexible with slow and fast internal motions (vibrations), the time step is as small as 0.5 fs (typically one twentieth of the highest frequencies of the vibrational modes). The cutoff radius is 10 Å and cubic periodic boundary conditions have been applied. Molecular Dynamics simulations have been performed in the NPT statistical ensemble (44) at constant standard pressure (1 bar) using a barostat relaxation time of 1 ps and a thermostat relaxation time of 0.2 ps. The temperatures range from 313 K to 493 K in steps of 20 K and the simulation lengths vary from 0.7 to 6 ns.

Results and Discussion

The intermediate scattering function S Q, t usually obtained from coherent neutron scattering experiments has been analyzed in order to investigate numerically collective dynamics.

The static structure factor S(Q) defined as S(Q, t = 0) is displayed in Fig. 5(a). A qualitative agreement with experimental data extracted from Ref. [START_REF] Champeney | [END_REF] is observed. This indicates that both real and model systems exhibit similar ordering. An interesting new point is the existence of a pre-peak (61) close to Q ≈ 0.6 Å-1 located below the first sharp diffraction peak at Q = Q max ≈ 1.4 Å-1 . This pre-peak seems to be also present in experimental data (65) at lower temperature. This feature has already been observed in some other hydrogen-bonded glass-forming molecular liquids such as metatoluidine and metacresol (66; 67). In such liquids, the pre-peak has been interpreted as a signature of structural heterogeneities corresponding to clusters of hydrogen-bonded molecules and forming a medium range order.

The dynamical structure factor S(Q, t) has been calculated from MD simulations at different temperatures and different Q wave-vectors (61) and is compared with the Neutron Spin Echo data (63) in Fig. 6 at Q=Q max . Despite slight discrepancies at the highest and the lowest temperatures, a reasonable agreement between experimental and numerical S(Q, t) is found. Upon cooling, S(Q, t) displays a decay usually observed in MD simulations (27) as well as in experiments on glass-forming molecular liquids (63; 68). The shape of S(Q, t) is similar to that drawn for glassy crystals in section 3. It is analyzed in terms of α-β relaxations in the framework of MCT (4). In the α long time relaxation regime, the characteristic time τ Q is defined as the time it takes for S(Q, t) to decay from 1 to 1/e. The characteristic time τ Qmax and the diffusion coefficient D obtained experimentally and numerically are shown in Fig. 7.

The characteristic time τ Qmax and the diffusion coefficient D are now consistent with previous experimental results (63; 62) and exhibit a non-Arrhenius dynamical behavior. According to MCT, any characteristic relaxation times should yield asymptotically the same temperature dependence ≈ (T -T c ) -γ . But an accurate determination of T c is difficult due to the limited temperature range investigated in Ref. (61). We can just check that the values (T c ≈ 288 K and γ ≈ 2.87) found in Ref. (69) Another intrinsic basis of MCT states that the behavior of any time-dependent correlators such as S(Q, t) describing the dynamics of the system is only controlled by its static density correlator S(Q). This statement has been nicely demonstrated both experimentally and numerically for the pure van der Waals glass-forming ortho-therphenyl system (70; 27). Fig. 5(b) shows the Q-dependence of the relaxation time τ Q , which is modulated in phase with S(Q). Moreover, the curve τ Q (Q) highlights the pre-peak (at ≈ 0.6 Å-1 ) and provides an estimation of the lifetime of the hydrogen-bonded molecules clusters, τ ≈ 80 ps.

Binary mixtures of Lennard-Jones atoms

From sections 3 and 4, we have seen that different molecular systems could exhibit similar dynamical behaviors upon cooling. However, the different kinds of interactions and the different sizes of the molecules make the comparisons ambiguous. Therefore, we have chosen to study the influence of the shape of the interaction potential on the dynamical properties of binary mixtures of LJ atoms.

Model

MD simulations have been carried out on binary particles mixtures with three different interaction potential (see Fig. 8). We have performed from 10 5 to 6 × 10 7 time steps, depending on temperature. All models are composed of 1500 uncharged particles (1200 of species A and 300 of species B). Particles interact via the (q, p) Mie potential having the form:

V (r) = E 0 (q -p) p r 0 r q -q r 0 r p . ( 4
)
The parameters r 0 and E 0 represent the position of the minimum of the well and its depth, respectively. The reduced LJ units (71) are used. The choice of q=12 and p=6 corresponds to the standard LJ potential (see Table 4) used by Kob & Andersen (K&A) , which has been studied extensively as a model glassformer (72; 73; 74; 75). For the purpose of investigating the change of dynamics with controlled change of V (r), we developed two other models by changing only the exponents, q and p, of the Mie potential for the A-A interaction. They are (q = 8, p = 5) and (q = 12, p = 11) and shown together with the (12,6) Lennard-Jones potential in Fig. 8. The (12,11), (12,6) and (8,[START_REF] Tölle | [END_REF] potentials are referred to as Model I, II and III, respectively. The well depth and the position of the minimum of V (r) are unchanged and the standard (12,6) Lennard-Jones potentials of the K&A model for the A-B and the B-B interactions are used in order to retain as much as possible the remarkable ability of the K&A model to form a glass upon cooling. The (12,11) model is more harmonic than the classical (12,6) potential, while the (8,5) potential is a flat well and exceedingly anharmonic.

Role of the shape of the interacting potential on the long time dynamics

From the self intermediate scattering function F S (Q 0 , t) at Q 0 = 2π/r 0 where r 0 is the position of the maximum of the first peak of the radial distribution function, we have determined a long relaxation time τ corresponding to the α process, when F S (Q 0 , t) has decayed by a factor of e -1 from the plateau-like region. In the MCT framework (see section 2), τ is proportional to (T -T C ) -γ . In order to check the validity of this scaling law, we have first determined the exponent b of Eq. 1 from the master-curve of F S (Q 0 , t) as displayed in Fig. 10. Another method would consist in using the factorization theorem in the β regime as described in (30). Knowing b, we have deduced a, leading to a value for the exponent γ. Once γ is known, we have fitted the α relaxation times τ (T ) -1/γ by a linear law of the form (T -T C ) in order to get the critical temperature T C of MCT. This linear behavior, only valid close to T C , is shown in Fig. 9 for the three studied systems. From the self intermediate scattering functions F S (Q 0 , t) and τ A (the α relaxation time τ of species A) at different temperatures for the three models, we can plot a master-curve obtained from the superposition of F S (Q 0 , t/τ A ) in the long time regime. The master-curve of Model III is displayed in Fig. 10. Upon lowering the temperature from the normal liquid state, we observe that the self intermediate scattering functions progressively collapse onto an apparent master curve. However, a detailed view of this thick curve, shown in the inset of Fig. 10, reveals that the different functions do not strictly superpose. Therefore, it is only approximatively possible to extract a temperature independent stretched exponent β using a Kohlrausch-Williams-Watt (KWW) fit of the master curve in a restricted temperature range. The same also stands for Models I and II. We have obtained β = 0.84 , 0.81 and 0.76 for Models I, II and III respectively. The error bar is about 0.1 due to the thickness of the master-curve.

5.3

Beyond MCT analysis, the binary LJ liquid-glass transition

5.3.1

Anharmonicity and intermolecular coupling. Since a strict superposition of F S (Q 0 , t/τ A ) functions does not hold, it is interesting to consider a temperature dependent stretched exponent β which is displayed in Fig. 11.

β is observed to continuously decrease from 1 at high temperatures to 0.69, 0.65 and 0.60 at very low temperatures below T C for Models I, II and III respectively. Moreover, β obtained from Model I is always found higher than β of Model II, itself higher than β of Model III. Therefore, the stretched exponent β decreases as anharmonicity increases. This behavior indicates that the anharmonicity of the intermolecular potential significantly controls the long range interactions between particles and therefore their collective relaxations. This result strongly supports the idea developed in the Coupling Model (76) that the intermolecular coupling measured by the parameter n = (1 -β) controls relaxations in liquids.

5.3.2

Link between short time and long time dynamics. We can focus now on the liquid-glass transition by looking at a characteristic parameter, that is the fragility m (defined below, see Eq. 5), and its sensitivity to the shape of the potential. The α relaxation time, τ , of all glass-forming liquids increases on cooling and becomes so long at some temperature T g that equilibrium cannot be maintained and the liquid is transformed to a glass as τ reaches commonly 10 2 s. This behavior is shared by very different glass-formers, but the temperature dependence of τ can differ greatly from one liquid to another (77; 78). Therefore, the change of log(τ ) with T g /T at T g /T =1, which is given by the steepness index or the fragility m, is defined by (79; 80; 81)

m = d log(τ ) d(T g /T ) Tg/T =1 . (5)
The values of m vary over a large range, from about 17 for strong glass-formers (like silica) to values as high as about 200 found for some glass-formers called fragile. Nevertheless, the microscopic origin of the large variation of the fragility m is not clearly understood yet. Several attempts have been made (see Ref. (82) for a list and references therein), one being the correlation between the fragility m and the vibrational properties of the glass at temperatures well below T g found recently by T. Scopigno et al. (83). The nonergodicity parameter, f (Q, T ) at T << T g is determined by vibrations and from inelastic X-ray scattering data, the temperature dependence of f 0 (T ) = f (Q → 0, T ) is well described by the following relation:

f -1 (Q → 0, T ) = 1 + α T T g . (6)
T. Scopigno et al. showed that m and α are proportional for many glass-formers. Before investigating the glassy state, we deal with the supercooled liquid state. In Fig. 12, log(τ A ) is plotted against T ref /T and the data of the three models show systematic change. It can be seen that the (2) gives values of m(τ A ) equal to 0.195, 0.241 and 0.405 for Models I, II and III, respectively. So, m(τ A ) increases with anharmonicity. We have also determined a steepness or fragility index, m(D A ) ≡ d log(1/DA) d(Tref /T ) as (T ref /T ) → 1 from the diffusion coefficient, D A , of species A calculated from the mean-square displacement < u 2 (t) > at long times when it reaches the linear t dependence. Again, m(D A ) increases monotonically in the order of Models I, II and III, or with anharmonicity.

So, the fragility m determined from different dynamical parameters and the intermolecular coupling n = (1 -β) increase with anharmonicity. We study now the very low temperature domain, well below T ref in order to concentrate on the vibrational properties of the three models in which relaxation is absent.

Therefore, the nonergodicity parameter f (Q, T ) determined from F (Q, t) is contributed entirely from vibrations. In all three models, the dependence of f 0 (T ) -1 on T /T ref shown in Fig. 13 is approximately linear, suggesting the following fitting law:

f 0 (T ) -1 = 1 + α T T ref , (7)
just similar to Eq. 6 used to represent the dependence of f 0 (T ) on T /T g of real glass-formers obtained by inelastic X-ray scattering (83). We see in Fig. 13 that the increase of f 0 (T )

), m(D A ), n = (1 -β))
and is even proportional to m(τ A ) (see the inset of Fig. 13) as shown in (83). At this point, we can state that the interacting potential is the origin of the correlation of α with m(τ A ) or m(D A ), and n = (1 -β), suggesting the same holds for real glass-formers. Fig. 14 summarizes the value of the fragility m estimated from different definitions for the three models identified by the intermolecular coupling parameter n = (1 -β). Whatever the definition used to calculate the fragility m, Fig. 14 shows that m of Model I is smaller than m of Model II, itself smaller than m of Model III, and that the increase of m is in relation with the increase of n. Moreover, before concluding, we observe the same trend for the fragility m determined by the relation presented in Eq. 8.

m = d log f dT T =T - ref (8)
It is then important to note that the parameter n, which is a measure of the intermolecular coupling due to the anharmonicity of the interaction potentials also controls the short time dynamics observed via the nonergodicity parameter f below T ref . This point is detailed in Ref. (82; 84). So, the shape of the interaction potential and particularly its anharmonicity is the key parameter, which controls the long time relaxations in supercooled liquids, the vibrations in glasses and which link them together.

Conclusion

Very different systems like orientationally disordered crystals, molecular and atomic glassforming liquids presented here exhibit similar dynamical behaviors on cooling. Other systems, not shown in the present manuscript, such as polymers or proteins exhibit the same behaviors. The main features discussed above are:

• the two step (α -β) relaxation which emerges at T A ,

• the existence of a critical temperature T C corresponding to an ergodic to non-ergodic transition,

• the fragility controlled by the anharmonicity (intermolecular coupling) of the interaction potential, • the link between short time and long time dynamics. The different systems considered above have been chosen to show the very diverse possibilities of dl poly.

The atomistic simulations are a very powerful tool in complement to theories and experiments, and permit us to find the microscopic mechanism responsible for macroscopic behaviors. A future challenge is a better description of the glassy state of glass-formers, which needs very long simulations that dl poly 3 could achieve.

F

  forming systems involving different degrees of freedom (translation and rotational): glassy crystals, glycerol solutions and binary mixtures of Lennard-Jones (LJ) atoms. Dynamical properties of these compounds have been particularly investigated in the ps-ns regime where a change in the nature of the dynamics (fast to slow dynamics) are expected to emerge in the [T A -T c ] temperature range.2 Two-step relaxation and Mode Coupling Theory

  , we have particularly shown from NMR and Raman experiments, and MD computer simulations that some predictions of the idealized version of the MCT (critical temperature T c and time scaling laws) with only TDOF were able to describe relatively well rotational dynamics of different plastic crystals.

3. 1 . 3

 13 Difluorotetrachloroethane. Difluorotetrachloroethane (DFTCE) is composed of simple molecules CFCl 2 -CFCl 2 close to dumbbells extensively used in MD calculations as prototype of molecular glassformers liquids

2

 2 Results and Discussion

  are compatible with the curves

  index m(τ A ) ≡ d log(τA) d(Tref /T ) as (T ref /T ) → 1, increases monotonically in the order of Models I, II and III. The values of m(τ A ) determined from this latter relation are 15.07, 18.57 and 26.58 for Models I, II and III. As a comparison, Sastry's method

Figure 1 .Figure 2 .Figure 3 .Figure 4 .Figure 5 .Figure 6 .

 123456 Figure 1. Reorientational correlation function C 2 (t) as function of time at different temperatures for a) CLADM, b) NBE and c) DFTCE. The dotted lines indicates the e -1 values.

Figure 7 .

 7 Figure 7. Temperature dependence of the relaxation time τ Q at Q=Qmax obtained from experiments (63) (full circle) and MD simulations (open circle). Temperature dependence of the diffusion coefficient D obtained from experiments (62) (full square), two MD simulations (open square (61), open triangle (85)).

Figure 8 .

 8 Figure 8. Potential V (r) determining the A-A interaction. The dashed curve is the(12,11) LJ potential for model I, the solid curve is the(12,6) LJ potential for model II and the dotted curve is the (8,5) LJ potential for model III.

Figure 9 .

 9 Figure 9. α relaxation time τ -1/γ of species A versus T /T C . Data can be fitted with a straight line ∼ (T /T C -1) for Model I ( ), Model II (•) and Model III ( ).

Figure 10 .Figure 11 .Figure 12 .Figure 13 .Figure 14 .

 1011121314 Figure 10. Check of the superposition principle of the self intermediate scattering functions F S (Q 0 , t/τ A ) for Model III. A detailed view is given in the inset.

  -1 with T /T ref is fastest for model III and slowest for model I, which means the slope α is largest for model III and smallest for model I. The parameter α increases with anharmonicity of the potential (like m(τ A

Table 1 .

 1 Parameters for the two-site CLADM model.

	o
	r
	P
	e e r
	R e v i e w
	O n
	l y

Table 2 .

 2 Simulation coefficients of NBE using a Buckingham potential φ(r) = -A/r 6 + B exp(-Cr).

	Site-Site A (kJ/mol Å6 ) B (kJ/mol) C ( Å-1 )
	C -C	2378.1	350142	3.60
	C -H	523.1	36701	3.67
	H -H	114.3	11112	3.74

Table 3 .

 3 Simulation coefficients of DFTCE using a Buckingham potential φ(r) = A exp(-ρ.r) -C/r 6 .

	Site-Site A (kJ/mol) ρ ( Å-1 ) C (kJ/mol Å6 )
	C -C	226307	0.288	2420
	C -F	196747	0.260	1168
	C -Cl	390940	0.284	3864
	F -F	171038	0.237	565
	F -Cl	320883	0.258	1808
	Cl -Cl	586389	0.284	5798

Table 4 .

 4 Parameters of the Lennard-Jones potentials that describe the three types of interactions in the K&A model(72; 73).

http://mc.manuscriptcentral.com/tandf/jenmol

Acknowledgments

The authors wish to acknowledge the use of the facilities of the IDRIS (Orsay, France) and the CRI (Villeneuve d'Ascq, France) where calculations were carried out. This work was supported by the INTERREG III (FEDER) program (Nord Pas de Calais/Kent).


Julie ; Diani
email: julie.diani@paris.ensam.fr

Bruno Fayolle

Lim

Lim Gilormini

Pierre ; Lim

P Gilormini

Study on the temperature dependence of the bulk modulus of polyisoprene by molecular dynamics simulations

Keywords: Molecular dynamics simulations, polyisoprene, bulk modulus, glass transition , compressibility Molecular dynamics simulations, Polyisoprene, Bulk modulus, Glass transition, Compressibility

HAL is

Introduction

The small-strain elastic mechanical behaviour of isotropic materials is fully characterized by two parameters. One can choose for instance, the Young's modulus and the bulk modulus. Dealing with amorphous polymers, a specific dependence of these two parameters vs. temperature is observed. The Young's modulus is known to drop by one to three orders of magnitude while temperature rises above the glass transition temperature [START_REF] Bueche | Physical properties of polymers[END_REF]. The bulk modulus is known to decrease by a factor of only 2 or 3 while temperature rises above the glass transition temperature [START_REF] Bicerano | Predictions of polymer properties[END_REF]. The temperature dependence of the Young's modulus has been explained well in the literature [START_REF] Treloar | Physics of rubber elasticity[END_REF]: above the glass transition, molecular mobility increases and enables entropic elasticity. However, the moderate drop of the bulk modulus is still not fully understood and the question of a possible entropy contribution to the bulk deformation is still open [START_REF] Matheson | Significance of entropic factor in mechanical deformation of polymeric glasses[END_REF]. Tabor's [START_REF] Tabor | The bulk modulus of rubber[END_REF] theoretical study of bulk modulus of rubber suggests that it arises from the Van der Waals interactions and that it is not related to the tensile modulus. Also, Theodorou and Suter's atomistic modelling [START_REF] Theodorou | Atomistic modeling of mechanical properties of polymeric glasses[END_REF] shows that entropy contributions to the elastic response to deformation can be [START_REF] Matheson | Significance of entropic factor in mechanical deformation of polymeric glasses[END_REF][START_REF] Matheson | Measurements of entropic factors in the mechanical deformation of polycarbonate and natural rubber[END_REF] reading of various experimental data of the literature seems to indicate that the entropy contribution cannot be generally neglected when considering the compressibility of polymeric glasses and rubbers. In this paper, we propose to address the question of the physical source of the bulk modulus of rubber (polyisoprene cis1-4) using molecular dynamics simulations.

F
Recently, the access to reliable affordable computers has contributed to the development of molecular dynamics simulations, and fully atomistic studies providing physical understandings on the mechanical behaviour of polymers are now possible. A number of studies have shown the ability of molecular dynamics to successfully reproduce the glass transition for amorphous polymers [START_REF] Rigby | Molecular dynamics simulations of polymer liquid and glass. I. glass transition[END_REF][START_REF] Han | Glass-transition temperatures of polymers from molecular-dynamics simulations[END_REF][START_REF] Yang | Molecular dynamics study of isobaric and isochoric glass transitions in a model amorphous polymer[END_REF][START_REF] Fukuda | Chain dynamics and conformation transition of cis-polyisoprene: Comparison between melt and subglasses state by molecular dynamics simulation[END_REF]. It is then reasonable to consider MD simulations to provide physical insight into the bulk modulus. This requires applying hydrostatic loadings, which are favourable for easy and fast MD simulations [12]. In such loadings, the number of degrees of freedom is reduced in comparison to uniaxial tension or shear. Using Materials Studio software [13], we built several virtual polyisoprenes and we applied hydrostatic compressions to the materials at various temperatures. During the simulations, the cell volume, pressure, and temperature were monitored, as well as the bonded, nonbond and kinetic energies. These data were used in a classical thermodynamics analysis to assess the physical source of the bulk modulus below and above the glass transition temperature.

General theory

We limit ourselves to hydrostatic loadings, which can result from an applied hydrostatic pressure p or from a volume contraction. The first law of thermodynamics relates a small change in internal energy U to small changes in entropy S and volume V, during quasi-static processes:

pdV TdS dU - = . ( 1
)
Considering the volume temperature T as two independent variables, Equation (1) leads to:

V V T S T T U ∂ ∂ = ∂ ∂ and T T V S T p V U ∂ ∂ + - = ∂ ∂ . ( 2
)
The isothermal bulk modulus is defined as:

T V p V V T B       ∂ ∂ - = ) , ( , (3)
and, by substituting Equation (2) in Equation (3) one gets

T T V S VT V U V V T B 2 2 2 2 ) , ( ∂ ∂ - ∂ ∂ = , ( 4
)
which allows to define the internal energy contribution and the entropy contribution to the bulk modulus as:

T U V U V V T B 2 2 ) , ( ∂ ∂ = and T S V S VT V T B 2 2 ) , ( ∂ ∂ - = (5)
In small strain elasticity, the bulk modulus refers to moderate volume contractions and is defined as

T B T V T B =
, where ) ( 0 T V corresponds to the volume of the material at temperature T and at atmospheric pressure. For various temperatures, we applied moderate volume contractions and used Equation [START_REF] Treloar | Physics of rubber elasticity[END_REF] to obtain the bulk modulus as a function of temperature. During MD simulations, the internal energy of the system is computed as the sum of the bonded energy, the nonbond energy and the kinetic energy:

c NB B K E E U + + = (6)
The nonbond energy consists of the Van der Waals energy, the electrostatic energy and the hydrogen bond terms, while the bonded energy accounts for the valence interactions, such as bond stretching, dihedral angle torsion, valence angle bending, inversion.... The force-field used in Material Studio is COMPASS [START_REF] Sun | COMPASS: An ab-initio force-field optimized for condensed-phase applications-Overview oith details on alkane and benzene compounds[END_REF], which is an abinitio-based force-field that has been validated on isolated molecules [START_REF] Yang | COMPASS force field for 14 inorganic molecules[END_REF][START_REF] Bunte | Molecular modeling of energetic materials: the parameterization and validation of nitrate esters in the COMPASS force field[END_REF]. In COMPASS, the Van der Waals interactions are defined by a 9-6 Lennard-Jones potential. Having access to the variation of the internal ernergy U as a function of V, we can assess the internal energy contribution

)) ( , ( 0 T V T B U
through Equation ( 5) and, finally, the entropy contributions to ) ( 0 T B is deduced.

Molecular dynamics simulations

Procedure

Fully atomistic simulations were performed on periodic systems of cis-polyisoprene. The virtual materials were built using an amorphous cell generator provided in Materials Studio, with periodic boundary conditions. Each parent cell contains 5 chains of 300 monomers. The cubic parent cells were built at a high temperature (above the expected glass transition) of 400 K to ensure a better stability. Each system was then submitted to consecutive, energy minimization, 1 ns-NVT MD simulation (number of atoms, volume and temperature remain constant in a NVT ensemble), and 1 ns-NPT MD simulations (number of atoms, pressure and temperature remain constant in a NPT ensemble). This provided reasonably stable amorphous structures at atmospheric pressure and 400 K. The structures were then cooled down to 100 K with a series of 25 K to 50 K jumps with intermediate NPT MD steps running for 500 ps with a time step of 1 fs. The duration of the simulations is constrained by the relatively high number of atoms (19 510) and by the fact that we use fully atomistic simulations (and not coarse grain simulations). Nevertheless, to ensure that 500 ps were long-enough duration, we monitored the cell specific volume during 1 ns for NPTs at a low and a high temperature. As shown in Figure 1, the volume has converged before 500 ps.

Six cells were used to calculate the glass transition of the virtual materials. From these six cells, two were submitted to volume contractions in order to study the Bulk modulus temperature dependence. For each temperature, both cells were submitted to moderate 1.5% and 3%-volume contractions in order to calculate ) ( 0 T B according to Equation (3). Cell 1 simulations were run using Berendsen [START_REF] Berendsen | Molecular-dynamics with coupling to an external bath[END_REF] barostat and thermostat, while cell 2 simulations used the Andersen [START_REF] Andersen | Molecular dynamics simulations at constant pressure and/or temperature[END_REF] barostat and thermostat. Both methods provided similar results. In addition, larger volume contractions were applied to both cells in order to better analyze the energy variations.

Results

During 1 bar-isobaric cooling, we observe a linear decrease of the specific volume at both high and low temperatures with decreasing temperature. The logarithm of the specific volume of the cells (volume per unit mass) is plotted vs. the temperature in Figure 2, so that the slope gives the coefficient of volumic thermal expansion directly. The change in the thermal contraction defines a clean glass transition at T g ~250 K, which is 50 K above experimental values reported in the literature (~203 K at 10 K/min) [START_REF] Tobolski | Polymer Science and Materials[END_REF] but is in good accordance with former MD studies on cis-polyisoprene [START_REF] Fukuda | Chain dynamics and conformation transition of cis-polyisoprene: Comparison between melt and subglasses state by molecular dynamics simulation[END_REF]. From Figure 2, the expansion coefficients above and below T g can be evaluated as In order to asses the material compressibility, we first submitted cell 1 to an increasing volume contraction at 350 K. The resulting pressure is plotted in Figure 3 as a function of the logarithm of the specific volume of the cell. We note that the material stiffens as the volume decreases, which indicates that compressibility decreases with volume contraction. This result is consistent with rubber behaviour reported in the literature [START_REF] Adams | The compressibility of rubber[END_REF]. As previously mentioned, the bulk modulus is related to the material response for small volume contractions and is defined by the linear approximation of p vs.

)) ( ln( T V according to Equation (3). Figure 3 illustrates how ) 350 ( 0 B is estimated to be 1460 MPa. For each temperature, the bulk modulus ) ( 0 T B is calculated using three states: atmospheric pressure and resulting volume ) ( 0 T V , 1.5%-volume contraction, 3%-volume contraction. In Figure 4, the values of the bulk modulus are plotted vs. temperature for cells 1 and 2, and a linear relation can be considered both above and below the glass transition temperature. Below T g , a rise in temperature produces a decrease of the bulk modulus twice as important as above T g . We obtained a value of 2100 MPa for ) 300 ( 0 B , which is in good agreement with the value of 2020 MPa reported in the literature [START_REF] Holownia | Effect of carbon black on Poisson's ratio of elastomers[END_REF]. Around the glass transition, ) ( 0 T B varies more significantly, and the ratio of the bulk modulus at 200 K (50 K below T g ) over the bulk modulus at 300 K (50 K above T g ) is close to 2. A change of such a magnitude through the glass transition is consistent with the experimental behaviour reported in amorphous polymers [START_REF] Bicerano | Predictions of polymer properties[END_REF]. Therefore, the bulk modulus temperature dependence presented by our virtual materials is representative of the behaviour of a real polyisoprene. Consequently, in an attempt to provide some physical understandings of the material compressibility, the system energy changes have been monitored during volume contractions at fixed temperature and during temperature changes at fixed volume. , where B k is the Boltzmann's constant, N is the number of atoms and T the temperature. Consequently, the kinetic energy remains constant during volume contractions at constant temperature. In Figure 5, the bonded energy is plotted vs. the cell specific volume for the three temperatures. In order to provide a simple read of the data in a single Figure for all temperatures, the bonded energy was normalized for each temperature by its maximum value, which was obtained for the highest contraction (or equivalently smallest volume). At 350 K, above T g , the bonded energy remains constant through volume contractions. At 150 K, the bonded energy fluctuates around a plateau value for moderate contractions, including those applied to determine the bulk modulus. Then, for severe volume contractions, the bonded energy increases drastically. At 250 K, the same trend for the bonded energy as at 150 K is recognized, though it is subtler: the bonded energy does not vary for moderate volume contractions and it increases evidently for high contractions. Therefore, it can be concluded that, at both high and low temperatures, the bonded energy does not vary significantly when volume contraction keeps within the range used to compute the bulk modulus.

As mentioned previously, three sources may contribute to the nonbonded energy: the Van der Waals interactions, the electrostatic forces and the hydrogen bonds. Since polyisoprene has no hydrogen bond and electrostatic forces are found to remain constant with temperature and volume contractions, the computed nonbond energy changes amount to Van der Waals energy changes only. The latter are shown in Figure 6 vs. the cell specific volume. For each temperature, the curve shape is reminiscent of a Lennard-Jones potential. The Van der Waals energy per unit mass varies with temperature, but Figure 7 shows that the three curves coincide reasonably if a suitable vertical translation is applied. This result allows the Van der Waals energy per unit mass to be written as:

) ( ) , ( ) , ( 0 0 T T V T E V T E VdW VdW - + = β , (7)
with 356 = β J.kg -1 .K -1 if 0 T is taken as 350 K. From Figures 5 and6, we conclude that when the virtual polyisoprenes are submitted to moderate volume contractions (which are the contractions that define the bulk modulus), only the Van der Waals energy varies significantly among the various terms contributing to the internal energy. Cell 2 was submitted to volume contractions at 400 K, 350 K and 150 K and lead to similar results than those presented above; they will be used to strengthen the analysis in the next section. The data obtained from the volume contractions at constant temperature will be useful to determine

)) ( , ( 0 T V T B U
and, by subtracting to ) ( 0 T B , to estimate the contribution of entropy to the bulk modulus.

Van der Waals contribution to the bulk modulus

Since the Van der Waals energy changes, in the range of volume contractions considered to compute the bulk modulus, leads the internal energy changes, Equation [START_REF] Tabor | The bulk modulus of rubber[END_REF] gives

6 T VdW T U U V E T V V U T V T V T B T B 2 2 0 2 2 0 0 0 ) ( ) ( )) ( , ( ) ( ∂ ∂ ≈ ∂ ∂ = = . ( 8
)
Figure 7 suggests that a simple equation should be able to approach the Van der Waals energy per unit mass.

Actually, the 9-6 Lennard-Jones potential used at the local scale to define the Van der Waals interactions suggests trying to apply the same function to the overall Van der Waals energy of the cell, and the Equation ( 7) would write:

) ( ) , ( 0 2 3 T T V c V a V T E VdW - + - = β . (9)
As shown Figure 8, Equation ( 9) does fit the Van der Walls energies for cells 1 and 2 relatively well when a=1.90 10 -3 J.m 9 . kg -4 , c=2.81 J.m 6 .kg -3 are used for cell 1 and a=2.0 10 -3 J.m 9 . kg -4 c=2.96 J.m 6 .kg -3 for cell2.

By substitution of Equation ( 9) in Equation ( 8), the following simple form is obtained for

) ( 0 T B U : ) ( 6 ) ( 12 ) ( 3 0 4 0 0 T V c T V a T B U - = (10)
For any temperature T, ) ( 0 T V can be evaluated from the two-part linear fit shown in Figure 2, and thus

) ( 0 T B U
can be computed from Equation ( 10) and compared to

) ( 0 T B
. This comparison is presented in Figure 9 for cells 1 and 2. For temperatures below the glass transition, it can be observed that the bulk modulus

) ( 0 T B is close to ) ( 0 T B U
, though the values of ) ( 0 T B are slightly underestimated in cell 2. Since the latter is computed from the Van der Waals energy only, which dominates the internal energy, no significant entropy contribution is to be noticed in the glassy state. In contrast,

) ( 0 T B U
overestimates the bulk modulus above the glass transition. This overestimation indicates that entropy does contribute to the bulk modulus above the glass transition. We note that the entropy affects the bulk modulus by reducing it (in the sense that

) ( 0 T B is lower than ) ( 0 T B U
). In Figure 9

Conclusion

Molecular dynamics simulations were performed on virtual polyisoprenes in order to investigate the physical source of the bulk modulus above and below the glass transition. The materials were submitted to moderate volume contractions at fixed temperature through NVT simulations, in order to estimate the bulk modulus dependence on temperature. The bulk modulus values thus obtained were shown to be in qualitative and quantitative good agreements with data reported in the literature. In particular, the bulk modulus showed a drop by a factor of about 2 at the glass transition. Investigating the internal energy change during volume contractions Waals energy by a 9-6 Lennard-Jones potential, and by applying a classical thermodynamic analysis, the internal energy contribution to the bulk modulus was estimated below and above the glass transition. This contribution was shown to reproduce fairly well the bulk modulus below the glass transition but to overestimate it above. In the latter temperature range, entropy changes lessen the bulk modulus by an amount that can be considered as temperature independent to first order.

  Study on the temperature dependence of the bulk modulus of polyisoprene by molecular dynamics simulations

  × 10 -4 K -1 , respectively, in qualitatively good agreement with values reported in the literature [19], = r α 6.16×10 -4 K -1 and = g α 2.07×10 -4 K -1 .

Cell 1

 1 was submitted to several volume contractions at three fixed temperatures 350 K, 250 K and 150 K and stabilized by NVT simulations. During these simulations, we have collected the changes of the three component energy (Equation (6)): the bonded energy, the nonbonded energy and the kinetic energy.Considering the latter, only a change of temperature affects K c since by definition

  may be considered as independent of temperature to the first order, with an average value of -400 MPa for cell 1 and -300 MPa for cell 2.

  temperature, we observed that for moderate contractions, sufficient to estimate the bulk modulus, the internal energy change arises from the Van der Waals energy only. Hence, by approaching the overall Van der

Figure 1 .

 1 Figure 1. Convergence of the specific volume of the virtual material vs. simulation duration during NPT ensemble.

Figure 2 .

 2 Figure 2. Estimate of virtual materials glass transitions and expansion coefficients through the plot of the logarithm of the cell specific volume vs. temperature.

Figure 3 .

 3 Figure 3. Resulting hydrostatic pressure vs. logarithm of the cell specific volume during volume contractions through NVT ensemble at 350 K.

Figure 4 .

 4 Figure 4. Cells 1 and 2 bulk modulus vs. temperature.

9 Figure 5 .

 95 Figure 5. Normalized bonded energy of cell vs. specific volume resulting from NVT simulations at constant temperatures 350 K, 250 K and 150 K.

Figure 6 .

 6 Figure 6. Van der Waals energy per unit mass of the cell vs. specific volume resulting from NVT simulations at constant temperatures 350 K, 250 K and 150 K.

Figure 7 .

 7 Figure 7. Superposition of the data plotted in Figure6showing the independence on temperature of the changes of the Van der Waals energy vs. volume changes.

Figure 8 .

 8 Figure 8. Fit of the cell Van der Waals energy variations with volume by a 9-6 Lennard-Jones potential: a)cell 1, b)cell 2.

Figure 9 .

 9 Figure 9. Comparison of the Van der Waals energy contribution to the bulk modulus with the material bulk modulus: a)cell 1, b)cell 2.

F 1 ) 1 )

 11 Cell specific volume (m 3 /kg) Van der Waals energy (J.kg -Cell specific volume (m 3 /kg) Van der Waals energy (J.kg -Page 17 of 18 http://mc.manuscriptcentral.com/tandf/jenmol

http://mc.manuscriptcentral.com/tandf/jenmol

Acknowledgments

This work was supported by grant JC05_43403 from the French Agence Nationale de la Recherche.