S Klotz

B Baptiste

T Hattori

S M Feng

Ch Jin

K Béneut

J M Guigner

I Estève

Synthesis and characterisation of a new graphitic C-S compound obtained by high pressure decomposition of CS 2

Keywords: graphite, high pressure, carbon disulphide, neutron scattering, Raman scattering

Carbon disulphide (CS2) is, together with its closest analogue CO2, one of the simplest molecular systems made of double covalent bonds. Under high pressure, the molecular structure is expected to break up to form extended crystalline or polymeric solids. Here we show that by compression at 300 K to ~10 GPa (100 kbar) using large-volume high pressure techniques, a sudden reaction leads to a mixture of pure sulphur and a well-defined compound with stoichiometry close to C2S which can be recovered to ambient pressure. We present neutron and x-ray diffraction as well as Raman data which show that this material

Introduction

The general trend in molecular systems under compression is that they solidify (in most cases crystallise) in the GPa range and at sufficient high pressure ultimately form extended crystalline or amorphous compounds as the molecules break up [START_REF] Hemley | Effect of Pressure on Molecules[END_REF]. These high pressure solids have vastly different electronic properties compared to the starting material and in some cases can be quenched to ambient pressure where they may have useful properties. A wellknown example is carbon-monoxide (CO) which transform under pressure of 5 GPa to a polymeric compound which can be recovered to atmospheric pressure and which is a highdensity energy material [2]. Another case is CS2, a colourless liquid with a strong odour which is isostructural to CO2. Following a reaction route indicated by Bridgman [START_REF] Bridgman | Freezing parameters and compressions of twenty-one substances to 50,000 kg/cm 2[END_REF], CS2 transforms at 3-5 GPa and ~250 °C into a black solid (« Bridgman black ») which can be recovered to ambient pressure and temperature. Analysis suggests that it consists of a disordered 3D network with the same stoichiometry as CS2 and highly crosslinked chains with both single and double C-S bonds [START_REF] Ochiai | Carbon dioxide and carbon disulphide as resources for functional polymers[END_REF][START_REF] Whalley | Structure of Bridgman's black carbon disulphide[END_REF][START_REF] Butcher | Formation and properties of the black form of carbon disulphide[END_REF][START_REF] Chan | Structure and electronic properties of solid polymeric carbon disulphide[END_REF][START_REF] Dias | Insulator-metal transition of highly compressed carbon disulphide[END_REF][START_REF] Bolduan | High pressure Raman study of solid CS2[END_REF][START_REF] Agnew | Pressure-and Temperature-Induced Chemistry of Carbon Disulfide[END_REF][11]. The material is a semiconductor [START_REF] Butcher | Formation and properties of the black form of carbon disulphide[END_REF] with a poor electrical conductivity of 10 -13 Ohm -1 cm -1 . Compression at 300 K (without heating) in a diamond anvil cell leads to a transition at 8-10 GPa which is governed by kinetics [START_REF] Bolduan | High pressure Raman study of solid CS2[END_REF], and it appears that the structure of the polymer depends on the specific experimental conditions it has formed [START_REF] Agnew | Pressure-and Temperature-Induced Chemistry of Carbon Disulfide[END_REF].

Here we report the transformation of CS2 in the 10 GPa range, at ambient temperature, using a large-volume ("Paris-Edinburgh") pressure device. This technique allows to load and compress ~100 mm 3 liquid samples, i.e. macroscopic quantities, between two opposed anvils and a metallic gaskets. It can be used for both in-situ characterisation by neutron and x-ray diffraction as well as ex-situ sample synthesis, see ref. [START_REF] Klotz | Techniques in High Pressure Neutron Scattering[END_REF] for a detailed description of the method. Our investigations involved two types of experiments: In a first measurement, the transformation of CS2 was observed and studied by in-situ neutron diffraction to 9 GPa. These revealed a sudden transition at ~ 8 GPa from a crystalline molecular solid to a disordered compound which can be recovered to ambient pressure. In a second type of experiments, the recoverability of the sample was exploited to produce in several loadings a considerable amount of material (~0.5 cm 3 ) at atmospheric pressure. The availability of such macroscopic sample quantities at ambient pressure, outside the pressure cell, then enabled us an in-depth characterisation of the reaction product by applying a variety of techniques, in particular x-ray and additional neutron diffraction, Raman scattering, infrared absorption, density and resistivity measurements.

Experimental

Neutron diffraction at ambient and high pressure

In-situ high pressure diffraction measurements were carried out at the PLANET beamline of the Japanese neutron Facility MLF at J-PARC, Tokai, Ibaraki [START_REF] Hattori | Design and performance of high-pressure PLANET beamline at pulsed neutron source at J-PARC[END_REF] using a VX4 Paris-Edinburgh press [START_REF] Klotz | Techniques in High Pressure Neutron Scattering[END_REF]. Approximately 30 mm 3 of CS2 purchased from Wako Pure Chemical Industries Co. Ltd (Japan) was loaded into double-toroidal sintered diamond anvils and null-scattering TiZr gaskets. The sample was then compressed at 300 K in steps of 1-2 GPa, up to ~8 GPa where it transformed suddenly. Pressures were determined to a precision of ± 1 GPa by the load on the anvils and a pre-determined calibration curve. Data in the crystalline phase (below 8 GPa) were analysed by Rietveld refinement using FullProf [START_REF] Rodriguez-Carvajal | Recent advances in magnetic structure determination by neutron powder diffraction[END_REF]. Once the transformation occurred, the sample was decompressed in steps to 0 GPa, again with data taken during the download.

Neutron diffraction patterns at ambient pressure from recovered and purified samples (see sections 2.2 and 2.3) were collected in a vanadium can. Pair correlation functions g(r) of the patterns S(Q) were obtained with the software package 0529_FourierTranslation.exe developed by Kyushu University/Dr. Kohara (SPring8) (Japan) using number densities of 0 =0.01419/Å 3 and 0 =0.0104/Å 3 for the un-annealed and annealed samples, respectively.

Ex-situ (off-line) sample production

Since the in-situ neutron data demonstrated that the sample obtained at high pressure can be recovered to ambient pressure, a series of samples were produced in off-line runs using singletoroidal tungsten carbide anvils [START_REF] Klotz | Techniques in High Pressure Neutron Scattering[END_REF] and copper beryllium gaskets, with a VX3 Paris-Edinburgh load-frame. This different anvil/gasket configuration allowed a much larger sample volume of initially 106 mm 3 , as well an easy extraction of the sample pellet from the metallic gasket.

Again, the transition occurred suddenly and was observable by and audible "click" and a drop in load at 72 ± 3 tonnes on the anvils, which corresponds to 7-8 GPa according to a calibration curve for such anvils. The transition is perfectly reproducible and was observed in at least 20 runs, half of them in runs concerned by this report. Samples were decompressed to ambient pressure at a speed of ~ 1 GPa/min. Liquid CS2 was purchased from ARCOS ORGANICS (99.9%) and Tokyo Chemical Industry (98%) and loaded into the cell using a previously described method [START_REF] Klotz | Techniques for neutron diffraction on solidified gases to 10 GPa and above: Applications to ND3 phase IV[END_REF]. The source of the samples had no observable influence on the outcome of the experiments.

Sample purification, determination of density and stoichiometry.

Visible observation of the sample after extraction from the pressure cell revealed an inhomogeneous compound of crystalline sulphur in a matrix of black material (see Results section). To remove the sulphur, the pellet was filed, then thoroughly grinded in liquid CS2 and stirred in a high-power ultrasonic bath to dissolve the sulphur. The turbid fluid was centrifuged at 7500 rpm to sediment the black micrometer-size particles, and the whole procedure was repeated a second time with fresh CS2. After drying in vacuum, the powder appeared homogeneous and no trace of crystalline sulphur could be detected neither visually nor by Raman scattering.

The average density of the inhomogeneous pellet was then measured by Archimedes' method using distilled water and an electronic balance model L432 from BEL Engineering (resolution:

1 mg) and its associated DENS-01 density kit. Given the initial mass and volume of the pellet, the mass of the powder after removal of sulphur as well as the known density of crystalline sulphur, the density and stoichiometry x and of the remaining homogenous material CxS could be easily determined to an accuracy of 2% and 5%, respectively.

Raman and infrared measurements

Raman experiments were carried out using a confocal in-house optical set-up using a Jobin-Yvon/Horiba HR-460 spectrometer equipped with a monochromator with 1500 grooves/mm and a Peltier-cooled Andor CCD. Raman scattering was excited by a continuous argon ion laser emitting at 514.5 nm and focused onto a spot of about 2 μm. Spectra were collected in backscattering geometry with maximal 7 mW laser power on the powder sample. Infrared measurements were performed on a Bruker IFS66V/S Fourier transform infrared spectrometer working in vacuum and aligned in transmission geometry. Spectra were acquired in the frequency region 600-8000 cm -1 with resolution of 4 cm -1 and accumulation of 10 scans using the mid infrared instrumental configuration with Globar, KBr and MCT as respectively source, beam splitter and detector. A homogeneous pellet of 13 mm diameter was prepared by mixing 1 mg powder sample with 200 mg grinded and dried KBr powder as transparent matrix material.

X-ray diffraction and gravimetric measurements

X-ray diffraction data of the purified samples were obtained on an Xcalibur-S 4-circle diffractometer from Rigaku (Agilent) with a Sapphire CCD detector and a Mo anode (17.4 keV), and the samples filled into a 1 mm kapton capillary. Diffraction patterns from the sample and the empty capillary were collected up to Q=15 Å -1 in transmission geometry and the azimuthal integration of the 2D patterns was done with the CrysalisPro software (CrysAlisPro 1.171.38.46, Rigaku Oxford Diffraction, 2015). After subtraction of the capillary contribution, the pair distribution function was obtained from the Fourier transform of the scattering function S(Q) using the PDFgetx3 program [START_REF] Juhás | PDFgetX3: a rapid and highly automatable program for processing powder diffraction data into total scattering pair distribution functions[END_REF]. Thermogravimetric measurements used a Q600 SDT from TA Instruments with Pt crucibles under nitrogen flow. The sample with initially a mass of approximately 10 mg was heated at a rate of 3 K/min from 20 °C to 700 °C.

Transmission and Scanning Electron Microscopy (TEM/SEM)

Transmission Electron Microscopy (TEM) measurements were carried out using a JEOL 2100F field emission gun instrument operating at 200 kV in combination with an UltraScan 4000 Gatan CCD camera and an energy dispersive x-ray detector with 140 eV resolution for elemental mapping. For Scanning Electron Microscopy (SEM) measurements we used a Zeiss Ultra 55 field emission gun SEM operated at 15 kV. Backscattered electron or secondary electrons modes were applied to investigate the morphology. Energy dispersive X-ray spectrometry (EDXS) point analysis was carried out at a working distance of 7.5 mm using an EDXS QUANTAX system equipped with a silicon drift detector XFlash 4010 (Bruker). Data were processed with the software Esprit (Bruker).

Electrical resistivity measurements

Four-probe DC resistivity measurements were carried out with a laboratory-made sample holder consisting of two electrically insulated Al-plates which compact the powder sample on a strip of 7 mm length and 0.3x1.0 mm 2 cross section. The sample holder was then attached and electrically contacted to a standard sample "puck" with its associated sample stick from Quantum Design. The measurements were carried out with a PPMS (Quantum Design) between 2 and 300 K using a current limit of 100 A and a cooling/heating rate of 1 K/min.

Varying the current showed ohmic behaviour of the resistivity.

Results

We first report the in-situ observation of the high pressure transformation of CS2 by neutron diffraction (Fig. 1). The crystalline phase below 8 GPa can be well fitted by the known Cmca structure [START_REF] Baenziger | Crystal Structure and Molecular Motion of Solid Carbon Disulfide[END_REF] with no detectable contaminant in the diffraction patterns. Rietveld refinements show an anisotropic compression of the structure with no change of the C=S bond length to within ± 0.01 Å (the precision of the technique) nor a change in the orientation of the molecule in the unit cell. At a load of 71 tonnes corresponding to 7.9 GPa, the load suddenly dropped by 1 tonne and the diffraction pattern revealed a transformation to a strongly disordered (amorphous) solid. Diffraction patterns were collected for the following 12 hours where no further change could be observed, indicating that the transition was complete. The sample was then decompressed in 4 steps (60, 40, 20, 5 tonne corresponding to an estimated 7.5, 6, 2 and 0 GPa, respectively) and no changes could be observed apart from a trivial shift of all features to larger d-spacings as the load (pressure) is decreased. The sample obtained at high pressure is hence recoverable to ambient pressure, with no apparent modifications except a decrease in density.

After these observations, the experiment was repeated with a different gasket-anvil setup allowing three times more volume per loading and an easy extraction of the recovered sample.

The aim was to produce larger sample quantities, at ambient pressure, which allow a detailed characterisation by various other methods which cannot be carried out in the pressure cell.

Runs using this setup gave reproducibly spherical black spheres of ca. 6 mm diameter inside the metallic gasket with no trace of remaining liquid (Fig. 2). Observation under a laboratory microscope shows an inhomogeneous sample with yellowish streaks within a matrix of black material. Micro-Raman measurements identified the yellow component as crystalline octasulphur, as expected, by comparison with the known Raman spectrum. Since we were interested only in the black homogeneous materiel, the sulphur was removed by dissolving the grinded powder in liquid CS2, as explained above. It was checked by Raman scattering that this procedure does not modify the remaining component. The density of the initial (inhomogeneous) pellet was determined by Archimedes' method (see above) and gave a value of 1.95 ± 0.03 g/cm 3 . From this and the mass of the powder after removing the crystalline sulphur, a composition of approximately C2S and a density of 1.76± 0.03 g/cm 3 was determined for the remaining compound. This value might be compared to the density of pure graphite (2.26 g/cm 3 ) and pure crystalline sulphur (2.07 g/cm 3 ). For simplicity, we will call this compound "C2S" in the following. The Raman spectrum (Fig. 3) of this purified black powder revealed immediately a nano-sized graphitic compound, i.e. a material made of sheets of carbon on a honeycomb lattice. Such material, which includes pyrolytic graphite and graphene, has been extremely well characterized by Raman scattering [START_REF] Nemanich | First-and second-order Raman scattering from finite-size crystals of graphite[END_REF][START_REF] Tuinstra | Raman Spectrum of Graphite[END_REF][START_REF] Gupta | Raman Scattering from High-Frequency Phonons in Supported n-Graphene Layer Films[END_REF]. Its characteristic fingerprints are: The G-peak at 1586 cm -1 , the D-band at 1350 cm -1 , the 2D peak at 2690 cm -1 as well as an additional feature at 2940 cm -1 which appears in material with large disorder. The ratio of the D and G peak varies with the particle size and can be used to determine the average linear dimension of the carbon sheets [START_REF] Tuinstra | Raman Spectrum of Graphite[END_REF] which gave ~25 Å in our case. In addition, the position of the G-band at The lower (blue) Raman spectrum are published data of disordered graphite [START_REF] Nemanich | First-and second-order Raman scattering from finite-size crystals of graphite[END_REF] (called "glassy carbon" in [START_REF] Nemanich | First-and second-order Raman scattering from finite-size crystals of graphite[END_REF]).

A significant difference compared to pure graphite or its various disordered forms is the existence of two intense peaks at 318 and 394 cm -1 which are absent in any of the latter materials. The peaks cannot be assigned to modes of remaining CS2, nor to vibrations of residual crystalline sulphur, nor to those of S8 or even S2 monomers which have no modes between 250 and 410 cm -1 [START_REF] Jackson | A universal chemical potential for sulfur vapours[END_REF]. Since SEM and TEM measurements (see further below) show the presence of only two elements, carbon and sulphur, these modes have to be ascribed to vibration of C-S bonds formed between the carbon sheet and sulphur atoms. In contrast to the Raman data, the measured infrared absorbance in the 500-7000 cm -1 range (inset fig. 3) is almost featureless showing only two weak, broad and distorted features of at 1200 cm -1 and 1575 cm -1 . The latter is reminiscent to one of the two IR active modes in graphite (868 cm -1 and 1588 cm -1 ) [START_REF] Nemanich | Infrared active vibrations of graphite[END_REF]. Note that single-layered graphene has no IR mode at all.

To invest the thermal stability of the compound, thermogravimetric measurements were carried out between 20 and 700 °C, see fig. 4. It is found that C2S starts to decompose at ~200 °C reaching the maximum speed of transformation at 260 °C. The transformation is complete at ~350 °C after which the sample has lost approximately 50% of its mass. Since the initial stoichiometry (C2S) as well as the change in mass is known, the annealed compound must still contain approximately 9 at% of sulphur, i.e. must have a composition of approximately C10S. From the measured net heat flow during the transition we estimate the total amount of energy absorbed by the process to 3.7 J, i.e. 30 kJ per mole of released sulphur. This number is typically in the range of a strong van der Waals or ionic binding, i.e. the released sulphur is certainly not covalently bonded to carbon. Next, Scanning and Transmission Electron Microscopy (SEM and TEM, respectively) were carried out which support the average chemical compositions of the samples cited above but reveal a large spread in stoichiometry: Fig. 5 shows representative SEM and TEM images which find that the sulphur composition in the un-annealed sample ("C2S") varies between typically 10 and 40 at% between different grains. In the annealed samples ("C10S") it is between approximately 5 and 20 at%, and this composition remained approximately constant even for samples heated to 700 °C. Clearly, about 10% of the initial 30 at% sulphur must be tightly bonded to carbon. It was also found that the un-annealed samples were prone to changes in their visual aspect during electron bombardment. This must hence be related to migration of loosely bonded sulphur, the same fraction which can be removed by heating.

The structure of the powder samples was then determined by neutron and x-ray diffraction applying a pdf analysis, as appropriate for disordered systems. We note that given a composition of initially C2S, the scattering of neutrons is almost completely dominated by carbon (the coherent neutron cross sections of C and S are 5.6 barn and 1.0 barn, respectively, i.e. S contributes less than 9% to the total intensity of the C2S sample), whereas for x-rays, sulphur (Z=16) governs the intensity compared to carbon (Z=6) (the scattering power for xrays is proportional to Z 2 ). The combined neutron and x-ray data have therefore the potential to pinpoint structural features in quite detail.

Figure 6 shows a neutron diffraction pattern of C2S and compares it to published neutron data of what was called "turbostratic graphite" [START_REF] Burian | Structural studies of multiwall carbon nanotubes by neutron diffraction[END_REF] (the terminology of disordered graphitic systems varies widely between authors). The agreement is remarkable and leaves no doubt on the 2D structure with its sp 2 carbon bonding. A pdf analysis (lower panel) confirms this conclusion: The peaks at 1.41, 2.45, and 2.83 Å, 3.74 Å and 4.24 Å can all be identified with the first four neighbour C-C distances of the 2D honeycomb lattice of graphite or graphene with 2.41 Å edge length. However, one notices the complete absence of a peak at 3.35 Å which would be expected for the inter-layer nearest neighbour C-C distance in graphite. In a diffraction pattern it corresponds to the most intense reflection, the 002. Fig. 6: Neutron diffraction pattern (a) compared to published neutron pattern of "turbostratic graphite" (b) [START_REF] Burian | Structural studies of multiwall carbon nanotubes by neutron diffraction[END_REF]. The g(r) obtained from pdf analysis of (a) is shown in (c), together with the geometry in graphite/graphene and its five nearest in-plane carbon-carbon distances.

Figure 7 shows corresponding x-ray patterns along with a similar pdf analysis. Both patterns are from the same sample as used in the neutron study, but one pattern was collected after annealing at 250 °C for 10 h in a vacuum furnace. The difference between the annealed and the non-annealed sample identifies therefore unambiguously distances relate to sulphur only.

The g(r) of the un-annealed sample shows a strong peak at 2.0 Å followed by a massif of several peaks between 3 and 5 Å. Other minor peaks, in particular at larger distances, may be artefacts of the data analysis given the limited Q-range. A comparison with the annealed sample reveals the disappearance of the most intense peak at 2.0 Å as well as a further peak at 3.3 Å. The remaining peaks are at 1.47 Å, 1.87 Å, 2.43 Å, 2.83 Å, 3.80 Å and 4.24 Å are all -with the exceptions of the 1.87 Å peak -to within ~0.05 Å the same as found in the neutron data analysis. This confirms hence once again a hexagonal graphitic lattice. Beyond any doubt, the two peaks at 2.0 Å and 3.30 Å in g(r) (arrows in Fig. 7c) must be related to sulphur, the fraction which is removed by heating to above 200 °C. The two distances are in fact exactly what one would expect for an S3 trimer with its well-known 107 degree S-S-S bond angle, a local geometry found in various sulphur allotropes [START_REF] Degtyereva | Vibrational dynamics and stability of the high-pressure chain and ring phases in S and Se[END_REF]. The peak at 1.87 Å must also be related to sulphur, but to the ~ 10at% fraction which is still present after annealing, and its value corresponds exactly to what is expected for typical covalent C-S bonds. It is this fraction of sulphur which must be responsible for the low-energy Raman modes at 300-400 cm -1 since these exist even in the annealed samples. Indeed, frequencies in this range are typical for C-S bending modes [START_REF] Sheppard | The vibrational spectra of some organic sulphur compounds and the characteristic frequencies of C-S linkages[END_REF][START_REF] Yuan | Far infrared study of carbon disulphide at high pressure[END_REF].

A closer look at the three sharpest features in the diffraction pattern of fig. 7a reveals that their positions agree well with those of "polymeric sulphur" [START_REF] Cataldo | A study on the structure and properties of polymeric sulfur[END_REF], in particular Tuinstra's "fibrous sulphur" labelled S [START_REF] Tuinstra | The structure of insoluble sulfur S[END_REF][START_REF] Stolz | The structural properties of liquid and quenched sulphur II[END_REF]These are highly disordered, non-equilibrium states, usually obtained by quenching sulphur from the liquid and which do not dissolve in carbon disulphide. But in contrast to such reported bulk polymeric sulphur we do not see the associated infrared signature [START_REF] Cataldo | A study on the structure and properties of polymeric sulfur[END_REF] in our samples, nor its decomposition [START_REF] Cataldo | A study on the structure and properties of polymeric sulfur[END_REF] at 120 °C. From this we conclude that the majority of sulphur present in our samples cannot be a separate phase coexisting in a mainly carbon-rich matrix.

We have also investigated the electronic properties of un-annealed (C2S) samples. Fig. 8 shows measurements of the dc-resistivity  and its temperature dependence between 5 and 300 K.

We find a conductivity (1/) of 0.1 Ohm -1 cm -1 at 300 K which is 13 orders magnitude higher than reported from Bridgman's black [START_REF] Butcher | Formation and properties of the black form of carbon disulphide[END_REF]. The compound is a semiconductor. From the slope of the ln() vs 1/T Arrhenius plot at high temperatures (above 200 K) we derive a gap of 46 meV (inset of Fig. 8). At lower temperatures the same presentation deviates significantly from being linear. However, the T-dependence can be well fitted over the entire temperature range by a variable hopping range model [START_REF] Mott | Conduction in non-crystalline materials[END_REF] given by the expression: = 0 exp(T0/T) p where p is a fitting parameter related to the dimension d of the hopping mechanisms by p=1/(1+d). We find p=0.43, i.e. d=1.3. This appears to be consistent with a structure where short range hopping occurs along a 1-dimensional C-C bond and where hopping over larger distances probe the 2-dimensional nature of the carbon sheets.

Discussion and conclusion

The emerging picture from these studies can hence be resumed as following: When liquid CS2 is compressed at 300 K it crystallizes at 1.3 GPa. Further compression to 7-8 GPa leads to a chemical reaction giving a mixture of pure crystalline sulphur and a compound with nanometric sp 2 -bonded carbon layers, approximate composition of C2S and a density of 1.76± 0.03 g/cm 3 , see Fig. 9 for an illustration. Most of the sulphur remaining therein is weakly (van der Waals) bonded to the carbon sheets, is locally in form of S3 trimers with some resemblance to a polymeric state and can be removed by heating to above 200 °C. The remaining 10 at% sulphur is strongly (covalently) bonded to the carbon backbone and gives rise to Raman modes in the 300-400 cm -1 range. It persists to temperatures of at least 700 °C. For chemical reasons, this type of sulphur is most likely bonded to carbon located at the edges of the carbon flakes. Annealing at T>200 °C removes the van der Waals bonded species, leaving only covalently bonded sulphur resulting in an average composition of C10S.

The question then arises what mechanism leads to the formation of CxS instead of "Bridgman black" as found in previous 300 K compression studies [START_REF] Dias | Insulator-metal transition of highly compressed carbon disulphide[END_REF][START_REF] Bolduan | High pressure Raman study of solid CS2[END_REF][START_REF] Agnew | Pressure-and Temperature-Induced Chemistry of Carbon Disulfide[END_REF][11]. For this we note that carbon disulphide has a positive formation enthalpy of +88.7 kJ/mole [START_REF] Naghavi | High-pressure layered structure of carbon disulphide[END_REF], i.e. the molecule is inherently metastable. This contrast with its closest analogue, CO2, which has a large negative formation enthalpy of -393 kJ/mol [START_REF] Naghavi | High-pressure layered structure of carbon disulphide[END_REF]. As a result of this, its polymerisation or partial decomposition is exothermic, and a reaction will self-propagate if the heat cannot be evacuated rapidly. This is likely the case in large-volume pressure apparatuses where the sample surface/volume ratio is small, in contrast to the geometry in diamond anvil cells [START_REF] Dias | Insulator-metal transition of highly compressed carbon disulphide[END_REF][START_REF] Bolduan | High pressure Raman study of solid CS2[END_REF][START_REF] Agnew | Pressure-and Temperature-Induced Chemistry of Carbon Disulfide[END_REF][11] where the anvils, in addition, are excellent heat sinks. The sample in our high pressure method is therefore likely to self-heat to beyond the stability of Bridgman black, potentially into the P/T region which currently is believed to be "carbon+sulphur" [START_REF] Agnew | Pressure-and Temperature-Induced Chemistry of Carbon Disulfide[END_REF] or "CS2+sulphur" [START_REF] Bolduan | High pressure Raman study of solid CS2[END_REF]. If this is correct, our work would mean that in the reaction product so far believed to be simply "carbon" is in fact a new carbon-sulphur compound with stoichiometry C2S.

Finally, this new compound may be compared to other known graphitic compounds. Given the absence of sharp Bragg peaks reminiscent to graphite, in particular the absence of a sharp 002 reflection (see above), the compound is not a canonical "graphite intercalated compound" (GIC) in the strict sense. However, CxS bears some resemblance with graphene oxide, though its structure is fundamentally simpler since it contains only carbon and sulphur with no oxygen functional groups such as hydroxyl. Its closest relative could be a compound termed "RGO/S" [START_REF] Zheng | In Situ Sulfur Reduction and Intercalation of Graphite Oxide for Li-S Battery Cathodes[END_REF], judged by the available x-ray, Raman and thermogravity data [START_REF] Zheng | In Situ Sulfur Reduction and Intercalation of Graphite Oxide for Li-S Battery Cathodes[END_REF]. This is a carbonsulphur composite obtained by reducing graphene oxide (GO) at 600 °C in sulphur vapour.

Since the compound appears to be a promising cathode material for Li-S batteries, our CxS material might as well be suitable for energy storage devices. Using the high pressure technique applied in this study, the material can be easily produced in cm 3 quantity and -apart from its academic interest -might have interesting other technological application similar to other nanometric 2D graphene-type materials [START_REF] Perrozzi | Graphene oxide: from fundamentals to applications[END_REF][START_REF] Smith | Synthesis, properties, and applications of graphene oxide/reduced graphene oxide and their nanocomposites[END_REF]. The availability of a material with high surface to weight ratio, with physical properties which can be tailored through its sulphur stoichiometry and which is easily dispersible in an aqueous alcohol solution could find applications for gas sensors, membranes and supercapacitors [START_REF] Perrozzi | Graphene oxide: from fundamentals to applications[END_REF][START_REF] Smith | Synthesis, properties, and applications of graphene oxide/reduced graphene oxide and their nanocomposites[END_REF].

Fig. 1

 1 Fig. 1 Neutron diffraction patterns of CS2 as a function of pressure. Patterns were shifted

Fig. 2 .

 2 Fig. 2. Left: Sample (black sphere in the centre) surrounded by the brownish metallic gasket,

Figure 3 :

 3 Figure 3: Raman spectra (main figure, upper red curve) and infrared absorbance (inset) of C2S.

Fig. 4 :

 4 Fig. 4: Thermogravimetric data showing the loss of mass (red line, left scale) and its derivative

Figure 5 :

 5 Figure 5: Representative TEM (a) and SEM (d,e) images of un-annealed (a,d) and annealed (e)

Fig. 7 :

 7 Fig. 7: X-ray diffraction pattern of (a) un-annealed and (b) annealed (250°C) samples. The

Figure 8 :

 8 Figure 8: Temperature dependence of resistivity of C2S. The main picture shows a fit to a

  Indeed, a flake with 25 Å diameter (see above) contains approximately 200 carbon atoms from which ~ 30 are on edge positions, i.e. 14%. For a flake with 40 Å diameter, approximately 10% are on edge positions. In conclusion, given the wide distribution in the size of the flakes, these numbers are consistent with the measured average stoichiometry of C10S.

Fig. 9 :

 9 Fig. 9: Structural model of CxS: Sulphur is attached to carbon sheets, either by van der Waals

Acknowledgments and Funding Information

S.K. is grateful to J. Biscaras (IMPMC) and D. Hrabovski (Sorbonne Université) for help in the conductivity measurements, S. Renaudineau (Sorbonne Université) for access to the thermogravimetric equipment, and F. Skouri-Panet (IMPMC) for using the centrifuge. The SEM facility at IMPMC was supported by Région Ile de France Grant SESAME 2006 NOI-07-593/R, Institut National des Sciences de l'Univers (INSU)-CNRS, Institut de Physique-CNRS, Sorbonne Université, and the French National Research Agency (ANR) grant ANR-07-BLAN-0124-01. The neutron scattering experiments were performed at the Japanese neutron spallation source MLF under proposals 2020B0188 and 2021I0011. S.K. wishes to thank Prof. J. A. Gonzalez Gomez (Universidad de Cantabria, Spain) for discussions on Raman data presented here.

Data Availability. All data are included in the manuscript and available from the corresponding author upon request and will be deposited on the CNRS cloud https://mycore.core-cloud.net/


Kateryna Goloviznina

Zheng Gong

Agilio A H Padua
email: agilio.padua@ens-lyon.fr

The CL&Pol polarizable force field for the simulation of ionic liquids and eutectic solvents

Molecular dynamics simulation, and theoretical chemistry in general, have been central tools in the study of ionic liquids and eutectic solvents, both of which are fluid phases dominated by ionic interactions and strong hydrogen bonds. These systems are still relatively recent, so the amount of experimental information given their enormous diversity is still scarce. Computational studies have produced several important discoveries, besides helping in the interpretation of experimental findings. Here we review the latest developments in the models describing the interactions and conformations of ionic liquid phases with atomic detail, in particular polarizable force fields. In these all-atom models, the response of electron clouds to their electrostatic environment is represented by induced dipoles, therefore moving beyond the pair-wise additivity of fixed-charge force fields. Most of the conceptual framework to develop polarizable force fields has been available and implemented in major molecular dynamics codes, so the main challenge to develop a useful polarizable model for ionic liquids and eutectic solvents is how to handle their diversity in chemical structures and interactions, in order 1 to build a model that can represent solutions, mixtures and interfaces with materials, and that is not too difficult to extend to new systems. As usually, force field development is not only science: it's an art! The ingredients of successful force fields are often the result of pragmatism and heuristics, rather than mathematical complexity and absolutely rigorous physics (although retaining a sound physical basis in the parameters definitely helps). We focus on the CL&Pol polarizable force field for ionic liquids, protic ionic liquids and eutectic solvents, improving on previous-generation fixed charge models. The CL&Pol force field is built as a detailed but extendable and transferable model offering more reliable predictions of thermodynamic, structural and transport properties, which should contribute to the advancement of the field and towards the design of better solvents, electrolytes, lubricants, etc.

Introduction

Room-temperature ionic liquids have been a particularly active field of research for the past 20 years, 1 in which molecular dynamics simulation was instrumental to predict and understand many physical and chemical properties. The atomistic force fields, or interaction potential models, are a fundamental ingredient of a molecular simulation study: if the underlying molecular model is not good, the results of the simulation will not be meaningful.

Presently, at the forefront of this topic, several groups are actively developing force fields that incorporate polarization effects explicitly, thus representing the response of electron clouds to local electrostatic environments. Even highly detailed force fields are still orders of magnitude faster in computational terms than quantum mechanics, allowing us to access longer size and time scales. In ionic media, electrostatic forces need to be described carefully and the new polarizable models are more advanced than the traditional, fixed-charge models still in widespread use to study aqueous phases, organic compounds, biomolecules, interfaces, and of course also ionic liquids.

Ionic liquids (ILs) are salts formed by large organic ions, with asymmetric, flexible molec-ular structures and delocalized electrostatic charge, all of these features contribute to hinder crystallization and thus lower the melting point to near room temperature. The archetypal ionic liquid is butylmethylimidazolium bis(trifluoromethanesulfonyl)amide,

[C 4 C 1 im][Ntf 2 ]
(the anion is also known as TFSI): the cation has its charge delocalized over an aromatic 5-member imidazolium ring and has one nonpolar alkyl side chain, making it asymmetric; the anion has the charge delocalized over the N and the four O atoms, and is highly flexible with many accessible conformations. The extent of the liquid range, and the thermal and chemical stability of this IL are remarkable. Ionic liquids in which the ions are an acid-base pair (with the equilibrium largely displaced towards the ionic species) form the class of protic ionic liquids (PILs), of which an important example is ethylammonium nitrate (EAN). Hydrogen bonding is a key aspect in PILs (and can be also in "aprotic" ILs). The availability of these ionic fluid phases near room temperature opens up vast fields of application, as solvents for reactions, separations or for the synthesis of materials, as electrolytes in energy-storage devices, as lubricants, etc.

When speaking of ionic liquids we also include deep eutectic solvents (DES), 2,3 a more recent class of environmentally acceptable solvents, which are mixtures of an organic salt with a molecular compound, typically with the role of H-bond donor. Several interesting organic salts that are abundant, low-cost and environmentally friendly, but solid, can be used to make liquids at room temperature taking advantage of the phase diagram of the mixture with a molecular compound, by staying close to the eutectic composition. The archetypal DES is a mixture of cholinium chloride with urea, both common and non-toxic compounds that are solid at room temperature, which yield a liquid mixture near the eutectic composition of roughly 1:2 molar ratio. The properties of the predominantly-ionic DES, like those of ILs, are determined to a large extent by Coulomb interactions and by H-bonds involving ions, with also contributions from non-polar groups. DES have many similar applications to those of ILs. From a molecular simulation standpoint, ILs 4-14 and DES [START_REF] Perkins | Molecular Dynamic Simulations and Vibrational Analysis of an Ionic Liquid Analogue[END_REF][START_REF] Perkins | Experimental and Computational Studies of Choline Chloride-Based Deep Eutectic Solvents[END_REF][START_REF] García | The impact of charges in force field parameterization for molecular dynamics simulations of deep eutectic solvents[END_REF][START_REF] Alcalde | An experimental and theoretical investigation of the physicochemical properties on choline chloride -Lactic acid based natural deep eutectic solvent (NADES)[END_REF][START_REF] Doherty | OPLS Force Field for Choline Chloride-Based Deep Eutectic Solvents[END_REF] can be modelled using the same framework and scales of description. Ionic liquids and DES are complex, structured phases, dominated by strong Coulomb interactions and with persistent medium-range ordering. As such, they pose several challenges to simulation methods and to the force fields, which have been highly successful tools to simulate aqueous solutions, organic liquids and biomolecules, and should in principle be equally applicable to ionic phases.

The slow dynamics (due to a relatively high viscosity) and long-range structural features of ionic liquids require simulations with consequent size and time scales (typically 6-10 nm and 10-50 ns), significantly larger than those needed for typical molecular liquids, including aqueous or organic phases. In order to attain longer scales at lesser computational cost, coarse-grained (CG) models become invaluable, allowing the study of systems with persistent long-range ordering or with slow dynamics, for example the mesophases of thermotropic ILs (ionic liquid crystals). [START_REF] Ji | Effect of the Chain Length on the Structure of Ionic Liquids: From Spatial Heterogeneity to Ionic Liquid Crystals[END_REF] Several coarse-grained force fields have been developed for ionic liquids, in which typically one, or a few, CG sites ("beads") bearing an electrostatic charge represent the ionic head-groups, and each CG site on the side chains represents a few CH 2 units. Most CG force fields are specific to certain families of ILs [START_REF] Roy | An Improved Four-Site Ionic Liquid Model[END_REF][START_REF] Fajardo | Molecular Dynamics Simulation of Imidazolium CnMIM-BF4 Ionic Liquids Using a Coarse Grained Force-Field[END_REF] and are not easy to combine with other models, but the Martini force field, [START_REF] Souza | Martini 3: A General Purpose Force Field for Coarse-Grained Molecular Dynamics[END_REF] which is fragment-based, has been recently extended to ILs, [START_REF] Vazquez-Salazar | Martini Coarse-Grained Models of Imidazolium-Based Ionic Liquids: From Nanostructural Organization to Liquid-Liquid Extraction[END_REF] providing a tool to simulate many different systems, since the Martini force field is designed to allow mixing of fragments.

Coarse-grained models retain the essential physics of a molten salt, namely the charge ordering an the ionic correlations. However, these models average-out many molecular degrees of freedom, hiding certain details of specific interactions, such as H-bonding sites, electrostatic charge distributions or conformational features. So-called united-atom models [START_REF] Liu | A Novel United-Atom Force Field for Imidazolium-Based Ionic Liquids[END_REF] only hide a few degrees of freedom corresponding to H atoms, which are fused into the heavy atoms, thus keeping many important features. In our view, keeping a high level of detail is important for many studies, namely to design ionic liquids as solvents, lubricants or electrolytes, through variations in their molecular structures that lead to subtle effects and finely-tuned properties. This is one of the aspects that make room-temperature ILs interesting in terms of their physical chemistry, and also for many applications.

On the smaller side of the multi-scale approach, in which the electronic structure is computed using quantum mechanics, only very small systems can be studied (a few dozen ion pairs) and during short trajectories. These ab initio molecular dynamics methods are very powerful to study local, quasi-chemical effects but do not reach the scales needed to evaluate energetic, structural and dynamic properties of the ionic liquid phases or their interfaces. Therefore, at present, all-atom molecular dynamics strikes the best compromise between computational cost, level of detail, and size and time scales attainable. [START_REF] Kirchner | Multiresolution Calculation of Ionic Liquids[END_REF] Since only very little experimental information was available on ionic liquids and DES, as this field evolved simulation has contributed to discover interesting features about these systems, as well as serving as a tool to interpret macroscopic properties from the molecular level. [START_REF] Kirchner | Multiresolution Calculation of Ionic Liquids[END_REF] Also, it makes sense to benefit from advanced computational methods to select fewer, and better targeted, experiments. For example, the nanoscale segregation between polar and non-polar domains in ILs, where medium-length alkyl side chains are segregated from charged head-groups due the strong cohesion between the latter, was seen in all-atom simulations [START_REF] Lopes | Nanostructural Organization in Ionic Liquids[END_REF] and later confirmed by X-ray scattering. [START_REF] Triolo | Nanoscale Segregation in Room Temperature Ionic Liquids[END_REF] This heterogeneous nature of ILs has important consequences for solvation. Our understanding the correlations in ionic motions, concerning how conductivity and diffusion are related, was also largely improved through MD simulation studies. [START_REF] Kashyap | How Is Charge Transport Different in Ionic Liquids and Electrolyte Solutions?[END_REF] Understanding ionic transport is key for the use of ILs as electrolytes.

The design of task-specific ionic liquids for CO 2 capture is another example of a research effort in which experiments and computational chemistry made progress together. [START_REF] Sheridan | Role of Molecular Modeling in the Development of CO2-Reactive Ionic Liquids[END_REF] The field of DES is newer, and simulating DES is in certain aspects more complex: they are by definition multicomponent mixtures containing both ions and neutral molecules, with strong H-bonds involving ions, which are challenging to model. When compared to ILs, the species in DES have more freedom to create solvation environments with local compositions that differ from the average composition. In purely ionic media this is more difficult because of the constraint of electroneutrality. In DES, simulations have shown local solvation environments that enhance the solubility of nonpolar, hydrophobic compounds. [START_REF] Hammond | Tuning the Solvation of Indigo in Aqueous Deep Eutectics[END_REF] Because of the enormous variety of chemical structures (and compositions) of ILs and DES, and in the logic of using computational methods to identify and design novel solvents or electrolytes, molecular force fields should aim for a wide applicability, both in terms of the functional groups they can represent and the properties that can be computed. Traditional fixed-charge force fields for organic compounds, such as OPLS-AA, [START_REF] Jorgensen | Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids[END_REF] AMBER/GAFF [START_REF] Wang | Development and Testing of a General Amber Force Field[END_REF] or CHARMM/CGenFF, [START_REF] Vanommeslaeghe | CHARMM General Force Field: A Force Field for Drug-like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Fields[END_REF] offer the advantage of vast parameter sets covering many functional groups, and also of a degree of transferability that allows them to be ported to new molecules or ions with just reasonable effort, following well-documented parameterization procedures.

Unfortunately, this class of fixed-charge models proved unable, in general, to represent both structural and dynamic properties of ILs correctly, likely because they were not developed for ionic phases. When using integer ionic charges, dynamics is systematically predicted to be too slow when compared to experiment, [START_REF] Tsuzuki | Molecular Dynamics Simulations of Ionic Liquids: Cation and Anion Dependence of Self-Diffusion Coefficients of Ions[END_REF][START_REF] Chaban | Polarizability versus mobility: atomistic force field for ionic liquids[END_REF] roughly by one order of magnitude in diffusion coefficient or viscosity. The heat of vaporization seems also to be somewhat overestimated.

Scaling-down the ionic charges by a factor of about 0.8 37 can improve the sluggish dynamics by reducing the ionic cohesive energy. Physical arguments in favor of scaling charges include accounting for charge transfer and for effective polarization effects, [START_REF] Maginn | Molecular Simulation of Ionic Liquids: Current Status and Future Opportunities[END_REF] the the basis of these arguments in bulk liquid phases are not totally sound. Unfortunately, with scaled charges densities tend to become too low, by up to 5 %, and in mixtures and solutions the balance between solute-solute, solute-solvent and solvent-solvent interactions may be modified, compromising the prediction of solubility, miscibility and phase behaviour. [START_REF] Cui | Influence of Charge Scaling on the Solvation Properties of Ionic Liquid Solutions[END_REF] Studying mixtures and solutions involves calculations for molecular or ionic species in different environments and a charge scaling factor optimized for a certain state may not be transferable to different environments. [START_REF] Maginn | Molecular Simulation of Ionic Liquids: Current Status and Future Opportunities[END_REF] Charge scaling also affects dielectric properties. [START_REF] Schröder | Comparing reduced partial charge models with polarizable simulations of ionic liquids[END_REF] These shortcomings have motivated the development of more sophisticated force fields, which include to some extent a representation of electronic degrees of freedom through explicit polarization effects, thus moving beyond pairwise-additive interaction models. This representation of electron clouds is, of course, hugely simplified when compared to electronic structure calculations (such as density functional theory) but the computational cost is kept of the order of traditional molecular dynamics (MD) with additive, fixed-charge models.

These developments for ILs have followed the evolution of polarizable force fields for organic compounds and biomolecules, namely the CHARMM family of models, [START_REF] Lemkul | An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications[END_REF] and they have been the subject of a recent review. [START_REF] Bedrov | Molecular Dynamics Simulations of Ionic Liquids and Electrolytes Using Polarizable Force Fields[END_REF] Here we will give an account of additional, newer developments in polarizable force fields that extend their applicability to protic ionic liquids and to DES, both the site of strong H-bonds involving densely-charged ions. It is curious that along the way, "old" tricks from simulation of high-temperature molten salts [START_REF] Wilson | Polarization effects in ionic systems from first principles[END_REF] have been brought in to solve some technical issues with the polarizable models. The developments reviewed here aim for a general, transferable model, suitable to handle the huge variety of ILs and DES, so that it can be used as a tool in molecular design of novel materials.

Representing Polarization Explicitly

Polarization (or induction) terms can be incorporated in MD force fields through several methods. Adding point induced dipoles [START_REF] Ren | Consistent treatment of inter-and intramolecular polarization in molecular mechanics calculations[END_REF][START_REF] Ren | Polarizable Atomic Multipole Water Model for Molecular Mechanics Simulation[END_REF][START_REF] Rick | Potentials and Algorithms for Incorporating Polarizability in Computer Simulations[END_REF] or Drude induced dipoles [START_REF] Lemkul | An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications[END_REF][START_REF] Drude | The theory of optics[END_REF][START_REF] Van Maaren | Molecular Dynamics Simulations of Water with Novel Shell-Model Potentials[END_REF][START_REF] Lamoureux | Modeling induced polarization with classical Drude oscillators: Theory and molecular dynamics simulation algorithm[END_REF][START_REF] Lamoureux | A simple polarizable model of water based on classical Drude oscillators[END_REF] to atomic sites are nearly equivalent methods, based on different ways to materialize the induced dipoles.

Point induced dipoles are just vector quantities assigned to the atomic sites, whereas Drude induced dipoles are formed by two particles of opposite charge connected by a spring. The response of the induced dipoles to the electrostatic environment is determined by the polarizability of the atomic sites, and so polarizable force fields have one additional type of parameter, but fortunately with a clear physical meaning. Only minor deviations between point and Drude induced dipoles have been documented, and that for highly polarizable atoms. [START_REF] Schmollngruber | Comparing induced point-dipoles and Drude oscillators[END_REF] Using fluctuating point charges on atomic sites [START_REF] Rappe | Charge Equilibration for Molecular Dynamics Simulations[END_REF][START_REF] Rick | Dynamical fluctuating charge force fields: Application to liquid water[END_REF][START_REF] Wu | The electronegativity equalization method fused with molecular mechanics: a fluctuating charge and flexible body potential function for [Emim][Gly] ionic liquids[END_REF][START_REF] Bauer | Recent Applications and Developments of Charge Equilibration Force Fields for Modeling Dynamical Charges in Classical Molecular Dynamics Simulations[END_REF][START_REF] Cavalcante | Polarizability effects on the structure and dynamics of ionic liquids[END_REF] provides an alternative method to represent polarization, with some drawbacks, such as inability to describe out-of-plane polarization and not being easy to mix polarizable and non-polarizable particles, and some advantages, namely that no additional particles are required. With fluctuating charges, electronegativity and hardness parameters are used to compute on-the-fly the atomic partial charges, which are not set a priori, and so charge-equilibration models can naturally render charge transfer between molecules. However, taming charge transfer requires additional mechanisms. [START_REF] Nistor | A Generalization of the Charge Equilibration Method for Nonmetallic Materials[END_REF][START_REF] Verstraelen | The Electronegativity Equalization Method and the Split Charge Equilibration Applied to Organic Systems: Parametrization, Validation, and Comparison[END_REF] Fluctuating-charge models have not been developed to cover as many molecular structures as induced-dipole models. The implementation of the Drude induced dipole method in widely used codes, including CHARMM, [START_REF] Lamoureux | Modeling induced polarization with classical Drude oscillators: Theory and molecular dynamics simulation algorithm[END_REF] NAMD, [START_REF] Jiang | High-performance scalable molecular dynamics simulations of a polarizable force field based on classical Drude oscillators in NAMD[END_REF] Gromacs, [START_REF] Lemkul | Implementation of extended Lagrangian dynamics in GROMACS for polarizable simulations using the classical Drude oscillator model[END_REF] LAMMPS [START_REF] Dequidt | Thermalized Drude Oscillators with the LAMMPS Molecular Dynamics Simulator[END_REF] and OpenMM, 62 makes this technique the most obvious choice for ionic liquids and DES.

Polarizable force fields have been proposed for some ionic liquids, [START_REF] Yan | Molecular Dynamics Simulation of Ionic Liquids: The Effect of Electronic Polarizability[END_REF][START_REF] Schröder | Simulating polarizable molecular ionic liquids with Drude oscillators[END_REF][START_REF] Nakano | A wave-function based approach for polarizable charge model: Systematic comparison of polarization effects on protic, aprotic, and ionic liquids[END_REF][START_REF] Salanne | From molten salts to room temperature ionic liquids: Simulation studies on chloroaluminate systems[END_REF] and almost none for DES, but only a few efforts to develop systematic and general polarizable models have been undertaken. [START_REF] Borodin | Polarizable force field development and molecular dynamics simulations of ionic liquids[END_REF][START_REF] Starovoytov | Development of AMOEBA Force Field for 1,3-Dimethylimidazolium Based Ionic Liquids[END_REF][START_REF] Vázquez-Montelongo | Current Status of AMOEBA-IL: A Multipolar/Polarizable Force Field for Ionic Liquids[END_REF][START_REF] Mcdaniel | Ab Initio Force Fields for Imidazolium-Based Ionic Liquids[END_REF][START_REF] Goloviznina | Polarizable Force Field for Ionic Liquids[END_REF] The APPLE&P force field by Borodin [START_REF] Borodin | Polarizable force field development and molecular dynamics simulations of ionic liquids[END_REF] was pioneering, with excellent predictions of thermodynamic, structural and transport properties of ILs. But extending and combining this force field to model other compounds or materials is difficult because of the uncommon functional form (exponential-6 potential and combining rules).

Besides, this is a proprietary model whose parameters are not available in the open literature. Among recent Drude-based models, the SAPT-FF force field [START_REF] Mcdaniel | Ab Initio Force Fields for Imidazolium-Based Ionic Liquids[END_REF]72 and the transferable CL&Pol force field [START_REF] Goloviznina | Polarizable Force Field for Ionic Liquids[END_REF][START_REF] Pádua | Resolving Dispersion and Induction Components for Polarisable Molecular Simulations of Ionic Liquids[END_REF][START_REF] Goloviznina | Extension of the CL&Pol Polarizable Force Field to Electrolytes, Protic Ionic Liquids, and Deep Eutectic Solvents[END_REF] incorporate information from accurate quantum calculations using symmetry-adapted perturbation theory (SAPT), a suitable quantum chemistry method to obtain interaction energies and their decomposition into repulsive, electrostatic, induction and dispersion terms. The SAPT-FF model [START_REF] Mcdaniel | Ab Initio Force Fields for Imidazolium-Based Ionic Liquids[END_REF]72 for ILs is developed from scratch with all intermolecular parameters fitted to SAPT energies. This model has been extended with a study of polyethylene oxide in an IL 75 and recently of urea, [START_REF] Jeong | A Transferable Polarizable Force Field for Urea Crystals and Aqueous Solutions[END_REF] which point towards an eventual extension to DES. The functional form of this force field consists of van der Waals interactions modelled by an exponential Born-Mayer repulsive term and a series of C n /r -n , n = 6, 8, 10, 12 powers 77 for the dispersion term. The case-by-case and costly SAPT parameterization, as well as the incompatibility of the functional form to that of widely used force fields are obstacles to transferability and a broad usage of this force field.

The CL&Pol model [START_REF] Goloviznina | Polarizable Force Field for Ionic Liquids[END_REF][START_REF] Pádua | Resolving Dispersion and Induction Components for Polarisable Molecular Simulations of Ionic Liquids[END_REF][START_REF] Goloviznina | Extension of the CL&Pol Polarizable Force Field to Electrolytes, Protic Ionic Liquids, and Deep Eutectic Solvents[END_REF] presented in this review is built according to a different strategy, favouring generality and transferability. It was developed from the existing fixed-charge CL&P model, [START_REF] Lopes | Modeling Ionic Liquids Using a Systematic All-Atom Force Field[END_REF][START_REF] Lopes | A generic and systematic force field for ionic liquids modeling[END_REF] one of the most widely used for ILs, that had been parameterized following the philosophy and functional form of OPLS-AA (and it is also compatible with AMBER/GAFF). As such, the CL&Pol model is relatively easy to extend, with applicability already going beyond ILs and DES, having been used to study for example alkalimetal electrolytes, dissolved gases, dyes and interfaces with nanomaterials. The qualities of the CL&Pol force field were demonstrated in the calculation of thermodynamic, structural and transport properties, yielding significant improvements when compared to the previous generation of fixed-charge models. The CL&Pol force field, including parameter database, examples and system-building tools for major codes, is available at public code repositories (https://github.com/paduagroup).

Building the CL&Pol Force Field

The CL&Pol polarizable force field is built on the basis of the CL&P fixed-charge model following several steps, illustrated in Fig. 1: 1. Adding Drude induced dipoles derived from atomic polarizabilities.

2. Scaling down the well-depth of the van der Waals interactions to avoid double-counting of induction effects (which were implicit in Lennard-Jones terms of the CL&P model).

This modification is performed based on SAPT calculations, but a computationally light scheme was developed to enable easy extension.

Slightly adjusting atomic diameters to improve the calculation of density.

The intramolecular bonded parameters and the atomic partial charges are taken from the CL&P values. The intramolecular terms describing covalent bond stretching and va-Figure 1: Layout of the procedure to obtain the CL&Pol force field from the original CL&P fixed-charge model. lence angle bending are taken whenever available from the extensive OPLS-AA [START_REF] Jorgensen | Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids[END_REF][START_REF] Watkins | Conformational Analysis and Liquid-State Properties from Ab Initio and Monte Carlo Calculations[END_REF] and AMBER/GAFF 79 parameter sets. If not available, the force constants and equilibrium geometries are calculated using the quantum chemical methods recommended for these force fields, thus maintaining compatibility with many molecular compounds. Thermodynamic and transport properties computed from condensed-phase MD, by averaging over long trajectories of many nanoseconds, are not sensitive to details of vibrational modes (whose periods are 6 orders of magnitude shorter). Certain ideal-gas properties would be sensitive to details of bonded terms. The torsion (or dihedral angle) potentials, which affect conformations, and the non-bonded Lennard-Jones (LJ) and Coulomb terms are, on the other hand, critical for condensed-phase simulations. So-called "proper" torsion potentials act between groups of 4 atoms, 1-2-3-4 bonded through 3 bonds, rotating around the central 2-3 bond, and are key to describe conformations. "Improper" torsion potentials are applied to groups of 4 atoms of which 3 of them are bonded to a central, planar atom, usually the 3rd, in order to stiffen out-of-plane bending of sp 2 carbon or nitrogen atoms, for example.

The dihedral angle potentials do not determine the conformations alone, since atoms belonging to the same molecule also interact through LJ and Coulomb potentials if they are sufficiently far apart (4 or more bonds), which represents steric hindrance effects and prevents overlaps. Therefore, an interdependence exists between the torsion potentials and the LJ and Coulomb parameters. The OPLS-AA procedure requires, first, that charges and LJ parameters be obtained, and only then dihedral angle coefficients be fitted to quantum mechanically calculated torsion energy profiles. The torsion parameters of both CL&P and CL&Pol are obtained in this consistent manner [START_REF] Lopes | Molecular force field for ionic liquids composed of triflate or bistriflylimide anions[END_REF][START_REF] Watkins | Conformational Analysis and Liquid-State Properties from Ab Initio and Monte Carlo Calculations[END_REF] in an attempt to describe ion conformations accurately.

Drude Induced Dipoles

A Drude induced dipole is formed by two point charges of opposite sign, a positively-charged Drude core (DC), at the center the atomic site, and a negatively-charged Drude particle (DP), connected to its core by a harmonic spring with an equilibrium distance of zero. The pair of charges gives rise to an induced dipole under an external electrostatic field, with magnitude determined by the Drude charges and by the force constant of the spring,

α = q 2 D k D (1)
We opt, as in the CHARMM polarizable force field, [START_REF] Lamoureux | Modeling induced polarization with classical Drude oscillators: Theory and molecular dynamics simulation algorithm[END_REF] to set the force constant of the DC-DP harmonic bond to k D = 4184 kJ mol -1 and obtain the Drude charge from the atomic polarizability. Atomic polarizabilities for ionic liquids derived from first principle calculations are available in the literature. [START_REF] Heid | Quantum mechanical determination of atomic polarizabilities of ionic liquids[END_REF][START_REF] Heid | Additive polarizabilities of halides in ionic liquids and organic solvents[END_REF] Because of the small atomic mass comparable to the mass of the DP and the low polarizability (low electron density), hydrogen atoms are treated as non-polarizable, with their polarizability merged onto the heavy atoms to which they are bound.

In a MD trajectory, the degrees of freedom (DOF) associated with Drude dipoles can be relaxed iteratively to their equilibrium positions at each time step, but this self-consistent relaxed-Drude method is slow. A faster alternative is to integrate the equations of motion of the DP by assigning a mass to them (usually 0.4 u, which is subtracted from the respective DC). The aim is for this extended Lagrangian method to generate a trajectory as close as possible to that of the relaxed Drudes, so in order to prevent kinetic energy draining into the DC-DP DOF, the relative motion of the DP with respect to their DC is thermostated at very low temperature, [START_REF] Lamoureux | Modeling induced polarization with classical Drude oscillators: Theory and molecular dynamics simulation algorithm[END_REF] typically of 1 K, with the remaining DOF being thermostated at the desired temperature. Therefore, the implementation of the Drude model in MD codes requires a dual-thermostat scheme, which has been implemented in several MD packages. [START_REF] Jiang | High-performance scalable molecular dynamics simulations of a polarizable force field based on classical Drude oscillators in NAMD[END_REF][START_REF] Lemkul | Implementation of extended Lagrangian dynamics in GROMACS for polarizable simulations using the classical Drude oscillator model[END_REF][START_REF] Dequidt | Thermalized Drude Oscillators with the LAMMPS Molecular Dynamics Simulator[END_REF][START_REF] Huang | Molecular dynamics simulations using the drude polarizable force field on GPUs with OpenMM: Implementation, validation, and benchmarks[END_REF] It was recently shown that heat flow to the cold Drude DOFs may lead to unbalanced temperatures between intramolecular and translational center-of-mass DOF. [START_REF] Son | Proper Thermal Equilibration of Simulations with Drude Polarizable Models: Temperature-Grouped Dual-Nosé-Hoover Thermostat[END_REF] To correct this bias, a triple temperature-grouped Nosé-Hoover thermostat has been implemented in the OpenMM [START_REF] Son | Proper Thermal Equilibration of Simulations with Drude Polarizable Models: Temperature-Grouped Dual-Nosé-Hoover Thermostat[END_REF] and LAMMPS [START_REF] Goloviznina | Extension of the CL&Pol Polarizable Force Field to Electrolytes, Protic Ionic Liquids, and Deep Eutectic Solvents[END_REF] codes. This improved thermostating scheme provides more accurate dynamic and structural quantities.

Preventing the Polarization Catastrophe

In a Drude polarizable force field such as CL&Pol the electrostatic potential between the induced dipoles includes a short-range damping function:

U elst (r ij ) = i j̸ =i q i q j r ij + i j T ij (r ij ) q i q D,j r ij + i j̸ =i S ij (r ij ) q D,i q D,j r ij (2)
where T ij and S ij are short-range damping functions, plotted in Fig. 2. These damping functions represent smearing of the point charges at short range, which has a physical sense [START_REF] Dequidt | Thermalized Drude Oscillators with the LAMMPS Molecular Dynamics Simulator[END_REF][START_REF] Van Duijnen | Molecular and Atomic Polarizabilities: Thole's Model Revisited[END_REF][START_REF] Noskov | Molecular dynamics study of hydration in ethanol-water mixtures using a polarizable force field[END_REF][START_REF] Wang | Development of Polarizable Models for Molecular Mechanical Calculations. 3. Polarizable Water Models Conforming to Thole Polarization Screening Schemes[END_REF][START_REF] Taylor | The effect of Thole functions on the simulation of ionic liquids with point induced dipoles at various densities[END_REF] when describing the interactions of electron clouds, and also help to prevent instabilities in the trajectories due to excessive correlations between nearby dipoles, or to rare events in which DP are captured by nearby positive particles, for example DC of other atoms, causing a "polarization catastrophe".

We use Thole damping functions, S ij , [START_REF] Thole | Molecular polarizabilities calculated with a modified dipole interaction[END_REF] (3), controlled by the atomic polarizabilities α and a universal parameter a ij = 2.6. 84

S ij (r ij ) = 1 -1 + s ij r ij 2 exp(-s ij r ij ), (3)
s ij = a ij (α ij ) 1/3 = (a i + a j )/2 [(α i α j ) 1/2 ] 1/3 (4)
For most aprotic ionic liquids, Thole damping is sufficient to avoid instabilities. But in strongly H-bonded systems, including protic ionic liquids and DES, and particularly in the presence of densely-charged ions such as alkali metal cations or halogens, the interactions between certain charges and induced dipole need also to be damped (or smeared) at short range. [START_REF] Goloviznina | Extension of the CL&Pol Polarizable Force Field to Electrolytes, Protic Ionic Liquids, and Deep Eutectic Solvents[END_REF] In the CL&Pol force field, this damping of charge-dipole interactions is used only to prevent instabilities of simulations. Inspired by the MD studies of molten salts [START_REF] Wilson | Polarization effects in ionic systems from first principles[END_REF][START_REF] Salanne | From molten salts to room temperature ionic liquids: Simulation studies on chloroaluminate systems[END_REF][START_REF] Jemmer | Dipole and quadrupole polarization in ionic systems: Ab initio studies[END_REF][START_REF] Salanne | Polarization effects in ionic solids and melts[END_REF][START_REF] Salanne | Including many-body effects in models for ionic liquids[END_REF] we used an adaptation of the Tang-Toennies 91 (TT) dispersion damping function,

T ij (r ij ) = 1 -c ij • exp (-b ij r ij ) 4 k=0 (b ij r ij ) k k! (5)
where the parameters b ij and c ij are set empirically to 4.5 and 1.0. [START_REF] Goloviznina | Extension of the CL&Pol Polarizable Force Field to Electrolytes, Protic Ionic Liquids, and Deep Eutectic Solvents[END_REF] The TT damping function serves as a safeguard against events causing trajectory instabilities, preventing DP from being captured by nearby positive atoms (Fig. 3) without affecting properties of systems in normal circumstances. The DP on the O atom of the anion leaves its DC and is pulled towards the H atom of the cation head group. Hydrogen atoms are not polarizable and so do not carry a DP.

Transferability by a fragment approach

The CL&Pol force field is designed to be extendable to new compounds with as much transferability of parameters as possible. Large ions are decomposed into smaller fragments for the purpose of deriving the force field parameters, and the chemical functional groups of the fragments are kept as much as possible unchanged across molecular or ionic entities. Typically, an ion is represented by a charged head group and one or several neutral side chains, for example, 1-butyl-3-methylimidazolium is split into 1-ethyl-3-methylimidazolium, over which the positive charge is distributed, and butane serves as template to the atoms in the side chain that are sufficiently far from the charged head-group. [START_REF] Pádua | Resolving Dispersion and Induction Components for Polarisable Molecular Simulations of Ionic Liquids[END_REF] We keep the ethyl group in the first fragment to emphasise that the imidazolium ring affects partial charges on atoms up to two bonds along the chain. [START_REF] Lopes | Modeling Ionic Liquids Using a Systematic All-Atom Force Field[END_REF] On the contrary, if the charge is not spread over a conjugated π-system, but localized on a particular atom with several side chains attached to it, such as in quaternary ammonium or phosphonium cations, a fragment with short methyl groups such as tetramethylammonium is used. Often anions, which are typically much smaller than cations in the most usual ILs, are kept as entire ions. Naturally, in cases such as alkylsul-fonatea or tosylate, long-chain carboxylates or phosphates, separate fragments are defined for the negative head group and for the neutral side chains or groups.

Once the general procedure is established, a balance between the number of fragments and their size should be kept, in order to avoid over-parameterization of the force field. Some cases merit a special treatment. Thus, in ethylammonium nitrate, an important protic ionic liquid, and in choline chloride, present in many DES, the cations are treated entirely and parameterized specifically, since they are still relatively small and this makes the force field more accurate for these compounds. In more general IL and DES families it is much more productive to use a transferable fragment approach, even if some specific details are less well reproduced.

Separating Polarization from Dispersion Interactions

The parameters describing van der Waals interactions in CL&Pol are taken from the nonpolarizable CL&P force field for ionic liquids, or from OPLS-AA for molecular compounds. In these fixed-charge models, polarization effects are already present, albeit described implicitly, lumped with dispersion terms in the non-bonded Lennard-Jones potential. Now in the CL&Pol force field, polarization (induction) effects are accounted for explicitly through the Drude induced dipoles. Therefore, the original CL&P LJ parameters need to be rescaled to remove double counting of polarization terms, so that in CL&Pol the LJ terms represents only London dispersion. We perform this through a scaling factor, k ij , which is applied to the well-depth, ϵ, of the LJ potential, corresponding to the fraction of dispersion contribution in the sum of dispersion and induction energies,

k ij = E disp E disp + E ind (6)
This scaling factor is evaluated per fragment, and functional groups present in the same ion can have different values, depending on the fragment they belong to. We aimed to have the best possible resolution between the different terms making up the intermolecular energies, specially the induction and dispersion components, so we chose SAPT, [START_REF] Jeziorski | Perturbation theory approach to intermolecular potential energy surfaces of van der Waals complexes[END_REF][START_REF] Hohenstein | Wavefunction methods for noncovalent interactions[END_REF] an expensive method which doesn't suffer from basis-set superposition errors, to evaluate dimer interaction energies. [START_REF] Goloviznina | Polarizable Force Field for Ionic Liquids[END_REF][START_REF] Pádua | Resolving Dispersion and Induction Components for Polarisable Molecular Simulations of Ionic Liquids[END_REF] Only dimers of fragments that can be found in first coordination shells, namely cation-anion, cation-neutral, anion-neutral and neutral-neutral, are considered. The nonbonded interactions between ions of the same charge are not modified due to the larger separations (second solvation shell) arising from electrostatic repulsion. The values of k ij are nearly 1.0 for apolar-apolar fragment dimers (e.g. two butane molecules representing alkyl chains) because the induction component is small compared to the dispersive one. When polar or charged fragments are involved, the values of k ij can be 0.4 or even less, meaning that the polarization terms are dominant compared to dispersion.

The high computational cost of SAPT calculations prompted us to devise a general scheme to predict the scaling factors on the basis of accessible molecular quantities, such as total charges, polarizabilities and dipole moments of the fragments. [START_REF] Goloviznina | Polarizable Force Field for Ionic Liquids[END_REF] Based on the classical equations of London dispersion and Debye induction, the scaling factor can be expressed as

k ij = 1 + c 0 r 2 ij Q 2 i α j + Q 2 j α i α i α j + c 1 µ 2 i α j + µ 2 j α i α i α j -1 (7)
where the α are the polarizabilities of the fragments, Q their net charges, µ their dipole moments. The distance r ij is the equilibrium separation of a given fragment pair, and bustness of the prediction is demonstrated by the reliability of our approach when extending it to new structures.

c 0 = 0.
Mixing CL&Pol with the other polarizable force fields, for example the SWM4-NDP water model [START_REF] Lamoureux | A polarizable model of water for molecular dynamics simulations of biomolecules[END_REF] or a polarizable CO 2 model, [START_REF] Becker | Polarizable Force Field for CO2 in M-MOF-74 Derived from Quantum Mechanics[END_REF] is rather straightforward. But in these cases, only the influence of the already polarizable molecules (water or CO 2 ) on the CL&Pol fragments should be taken into account to scale the LJ terms (coming from CL&P),

k ij = 1 + c 0 r 2 ij Q 2 i α i + c 1 µ 2 i α i -1 (8)
where the i index corresponds to the already polarizable molecule. In the case of a water molecule, for example, this expression depends only on its dipole moment and molecular polarizability, [START_REF] Lamoureux | A polarizable model of water for molecular dynamics simulations of biomolecules[END_REF] resulting in a value of 0.72 for any dimer formed with SWM4-NDP water.

The procedure described so far to generate the CL&Pol parameters, by adding Drude induced dipoles and scaling down the well-depths of the LJ potentials, was seen to lead to systematically slightly lower densities by a few percent on the training set of ionic liquids. [START_REF] Goloviznina | Polarizable Force Field for Ionic Liquids[END_REF] This does not seem to be a spurious effect of the dual thermostat used to handle the Drude DOF, as has been reported in the literature, [START_REF] Son | Proper Thermal Equilibration of Simulations with Drude Polarizable Models: Temperature-Grouped Dual-Nosé-Hoover Thermostat[END_REF] since for the CL&Pol force field both classical and temperature-grouped Nosé-Hoover thermostat give comparable density deviations [START_REF] Goloviznina | Extension of the CL&Pol Polarizable Force Field to Electrolytes, Protic Ionic Liquids, and Deep Eutectic Solvents[END_REF] as seen in Fig. 5. Therefore, we interpret the lower densities as due to the procedure to parameterize the CL&Pol model, and opted to nudge the values of the LJ diameters σ by a universal factor of 0.985, i.e. a 1.5 % correction, in order to compensate for the density deviations. [START_REF] Goloviznina | Polarizable Force Field for Ionic Liquids[END_REF] For some particular cases, such as ammonium and phosphonium cations with several long alkyl side chains, for example trihexyltetradecylphosphonium (P 6,6,6,14 +), the correction needs to be more consequent, with σ scaled by 0.950 in order to keep density deviations within -2.5-3.0 %. This seems to be an issue arising from the very long side chains rather than from the ionic or polar fragments, and likely could benefit from sourcing the alkyl chain parameters from a native polarizable force field for alkanes.

-4

-2 0 2 4 [C 2 mim][DCA] [C 4 mim][BF 4 ] [C 4 mim][DCA] [C 4 mim][NTf 2 ] [C 6 mim][DCA] [SiCSiCmim][DCA] [SiCSiCmim][OTf] [SiCSiCmim][NTf 2 ] [C 8 mim][NTf 2 ] [Me 4 C 5 mim][NTf 2 ]
NH tgNH density overestimation (3.5 %) and "freezing" of the system, meaning the strong H-bonds involving densely-charged ions were not modelled correctly. The absence of a repulsive core on a hydrogen atom was allowing the positive charge of the H to approach Cl -and its DP, even when using the TT damping function. Note that H atoms bonded to O or to N in the OPLS-AA force field do not carry their own LJ sites; they are simply point charges embedded in the LJ site of the corresponding heavy atom (this is also the case in many water models).

ρ dev / % s σ = 1.000 s σ = 0.985
During the development of CL&Pol, attempts including a small LJ site on H atoms were made, to prevent H-bonds from becoming too short. But a satisfactory compromise between preventing instability in the Drude dipoles and maintaining correct structural features was not found. This led us to revise the source of the CL&P chloride parameters, which had been obtained from the crystal structures for the few available ILs more than 15 years ago. [START_REF] Goloviznina | Extension of the CL&Pol Polarizable Force Field to Electrolytes, Protic Ionic Liquids, and Deep Eutectic Solvents[END_REF] This leads to density, viscosity, diffusion coefficients, and also X-ray and neutron structure factors from polarizable MD simulations that are in good agreement with experimental data. [START_REF] Goloviznina | Extension of the CL&Pol Polarizable Force Field to Electrolytes, Protic Ionic Liquids, and Deep Eutectic Solvents[END_REF] It also produces linear H-bonds, supporting this hypothesis about the structure of EAN. Given the nature of the electrodes and of the IL electrolytes, it would be interesting to treat both using polarizable models, to representing more faithfully the interactions, ordering and dynamics in the interfacial layers, that are not yet fully understood. 111 But in most published studies using simulation, the two partners, electrode material or ionic liquid, are not treated with equivalent levels of detail, thus some important effects may be missed. We think in the future, with the availability of polarizable force fields such as CL&Pol, the community working on electrochemical devices will opt to model all-atom, fully polarizable systems. The CL&Pol force field has been used to study substituted imidazolium dicyanamide ILs at charged interfaces of molybdenum disulfide, MoS 2 , a two-dimensional material of the family of transition metal dichalcogenides. MoS 2 is a 2D semiconductor, interesting to make ionic liquid-gated devices, such as transistors in which an electrolyte is used to polarize the material and to modulate the density of charge carriers. 112,113 A parallel-plate capacitor was setup, as shown in Fig. 7, with the ionic liquids as the electrolytes and the image-charge constant-potential method representing the polarization of the MoS 2 surfaces, at imposed electrostatic potential. The effect of the side-chain length on the cations and of different substituents, namely aromatic and hydroxide groups, on the structure of the electrical double layer was analysed, providing details and insights about the orientation of ions and also on the capacitance and on the dynamics of charging the interface.

Conclusions

Polarizable force fields offer a productive and insightful framework to study systems with ionic liquid phases, be it as solvents, electrolytes, lubricants, confined liquids or coatings, among other possible roles. At present all-atom molecular dynamics simulation strikes an interesting compromise between level of detail and the system size and time scales attainable with reasonable computational cost. This is the mid-level, between quantum chemistry at the bottom, useful to investigate chemical or quasi-chemical phenomena (reactivity, charge transfer) at a local molecular scale, and coarse-grained models at the top, which are suitable for slow, long-range ordered systems such as ionic liquid crystals or polyelectrolytes. We learned over almost two decades that ionic liquids are sophisticated, their balance between Coulomb and van der Waals forces resulting in non-trivial structural and dynamical features.

Eutectic solvents, mixtures of salts with molecular compounds, are also complex with strong hydrogen bonds playing a major role (as they do in protic ionic liquids). The availability of the general force fields such as the CL&Pol model, built to be extendable to new compounds and systems without the need for case-by-case parameterization, is expected to contribute to advances in understanding the the physical chemistry of ionic fluids and their interfaces, thus enabling their rational design as better solvents or electrolytes.

Figure 2 :

 2 Figure 2: Short-range damping of Drude dipole-dipole interactions using a Thole function, S ij (a = 2.6, α = 1.5 Å 3 ) and of charge-Drude dipole interactions using Tang-Toennies functions, T ij (b = 4.5 Å -1 , c = 1.0).

Figure 3 :

 3 Figure 3: Illustration of the polarization catastrophe in ethylammonium nitrate, a protic IL.The DP on the O atom of the anion leaves its DC and is pulled towards the H atom of the cation head group. Hydrogen atoms are not polarizable and so do not carry a DP.

  25 and c 1 = 0.11 are adjusted coefficients to a set of SAPT-derived k ij values covering charged and neutral fragments, 71 as shown in Fig. 4. The root-mean square error does not exceed 0.08, which we find quite good, since conformations and local geometry minima also cause uncertainties in the SAPT values. Recently, the scheme was applied to dimers composed of new tetramethylphosphonium (P + 1111 ) and cholinium (Ch + ) cations, tetracyanoborate (TCB -) and 2,2,2-trifluoro-N-(trifluoromethylsulfony (TSAC -) anions, dimethyl ether (Me 2 O) and tetramethylsilane (TMS) molecules, and the ro-

Figure 4 :

 4 Figure 4: Predicted and SAPT-calculated k ij scaling factors. Empty points are those used in the original CL&Pol development, 71 whereas filled points are more recent involving new fragments.

Figure 5 :

 5 Figure 5: Predicting the density of ILs before and after LJ σ correction by 0.985. NH denotes the dual Nosé-Hoover thermostat, tgNH the temperature-grouped Nosé-Hoover thermostat; SiCSiC corresponds to dimethyl((trimethylsilyl)methyl)silyl)methyl, Me 4 C 5 to 2,2,4,4-tetramethylpentyl side chains in imidazolium cations.

Figure 6 :

 6 Figure 6: Two hydrogen bonding patterns in ethylammonium nitrate. A linear N H • • • O bond obtained from QM calculations in condensed phase 97,101,102 and the CL&Pol force field 74 (left) and a bent, shared N H • • • O bond obtained with the empirical potential refinement method from X-ray and neutron scattering 98-100 (right).

Figure 7 :

 7 Figure 7: Constant-potential simulation of a substituted imidazolium dicyanamide ionic liquid electrolyte near charged MoS 2 surfaces, with the polarization of the material represented by the image charge method.

  [START_REF] Bodo | A prototypical ionic liquid explored by ab initio molecular dynamics and Raman spectroscopy[END_REF] As said above, EAN is such an important liquid due to its network of hydrogen bonds, and this status grants a specific treatment in force field development (much like for water). Coincidentally, there have been contradicting analyses in the literature about the nature of the H-bonds in EAN. In one view, based on the empirical potential refinement method from X-ray and neutron structure factors, a number of authors[START_REF] Hayes | The Nature of Hydrogen Bonding in Protic Ionic Liquids[END_REF][START_REF] Song | Structural Heterogeneity and Unique Distorted Hydrogen Bonding in Primary Ammonium Nitrate Ionic Liquids Studied by High-Energy X-Ray Diffraction Experiments and MD Simulations[END_REF][100] postulated bent N H • • • O bonds, with an O atom from nitrate shared between several ammonium H. A different view, of linear H-bonds involving one ammonium H was advanced by other groups, 97,101,102 from DFT trajectories in condensed phase. The two patterns are shown in Fig. 6. The interpretation of the H-bond geometry from structural data seems to be highly sensitive to the underlying atomistic model, 103 so this is an open question. Going back to the extension of CL&Pol

	7,9 At
	the time this was thought to better represent the fully ionic nature of ILs than parameters
	for ions derived from hydration properties. But for strongly H-bonded systems, chloride
	parameters from OPLS-AA, which were obtained from aqueous solutions, 96 much improve
	the situation. Indeed in OPLS-AA Cl -has a larger LJ diameter (4.41 Å instead of 3.65 Å
	in the initial CL&P). Adopting OPLS-AA as a source of parameters for halides leads to a
	good prediction of equilibrium and transport properties of DES, such as cholinium chloride
	+ urea, cholinium chloride + ethyene glycol and tetrabutylammonium bromide + ethylene
	glycol, providing us with a suitable polarizable force field for DES.
	Quite similarly, in protic ionic liquids such as ethylammonium nitrate (EAN), a similar
	issue related to the shortening of N H • • • O distance was seen. There are three "naked-

charge" H atoms in the ammonium head group and their interactions with a nitrate anion led to an extremely short N H • • • O distance of 1.2 Å, comparable in length to an intramolecular N H bond, whereas ab initio calculations give a N H • • • O distance of 1.9 Å. to EAN, correct H-bond distances are obtained by increasing σ NO slightly.

Acknowledgement

The work presented here was supported by the LIQUI2D research project from the French ANR agency (ANR-18-CE09-0018) and by an IDEX Lyon Fellowship (ANR-16-IDEX-005).