Marek Szczerba

Andrey G Kalinichev

Intercalation of ethylene glycol in smectites: Several molecular simulation models verified by X-ray diffraction data

Keywords: Ethylene glycol, Smectite, Molecular dynamics, X-ray diffraction

Organo-clays represent a special challenge for molecular simulations, because they require accurate representation of the clay part and the organic/aqueous part of the model system and an accurate representation of their interactions. Due to a broad diversity of available force field models, it is an important question which sets of parameters will best suit the molecular modeling of the organo-intercalated smectites. To address this question, we use the structure of ethylene glycol (EG)-smectite complex as a testing model, because the intercalation of EG in smectites provides a stable interlayer complex with relatively constant basal spacing.

Three smectite samples with substantially different layer charge and charge localization were selected for X-ray diffraction (XRD) measurements. Their molecular models were built and molecular dynamics simulations were performed using various combinations of the organic force fields (CGenFF, GAFF, CVFF and OPLS-aa) with ClayFF and INTERFACE force fields used to describe smectites. The simulations covered a range of different EG and water contents in the selected smectites. For every structure, the density distribution of interlayer species along the direction perpendicular to the layer plane was calculated and then used to optimize the X-ray diffraction patterns for these simulated models.

A comparison of these results with experimental XRD patterns shows very large discrepancies in the obtained structures and basal spacings for different layer charges as well as for different force fields and their combinations. The most important factor affecting the accuracy of the calculated diffractograms is found to be the selection of the clay mineral force field parameters. The second important conclusion is that a slight modification of the basal oxygen parameters for non-electrostatic interactions (increase of their effective atomic diameters) may be a simple and straightforward way to greatly improve the agreement of the modeled diffractograms with experiments, especially for high-charge smectites. Generally, among the organic force fields the least accurate results are obtained with CGenFF. For unmodified ClayFF, its combination with GAFF gives the best results, while the two other sets (OPLS-aa and CVFF) give the best results in combination with ClayFFmod. The INTERFACE and INTERFACEmod force fields are found to be producing much better results for low-charge montmorillonite than for high-charge smectites.

INTRODUCTION

Recent growing interest in organo-clay nanocomposites is clearly reflected in the development of computational molecular modeling approaches to quantitative understanding of the structural and dynamic behavior of such systems (e.g. [START_REF] Zeng | Molecular dynamics simulation of organic-inorganic nanocomposites: Layering behavior and interlayer structure of organoclays[END_REF][START_REF] Suter | Computer simulation study of the materials properties of intercalated and exfoliated poly (ethylene) glycol clay nanocomposites[END_REF][START_REF] Suter | Rule based design of clay-swelling inhibitors[END_REF][START_REF] Suter | Chemically specific multiscale modeling of clay-polymer nanocomposites reveals intercalation dynamics, tactoid self-assembly and emergent materials properties[END_REF][START_REF] Greathouse | Interaction of Natural Organic Matter with Layered Minerals: Recent Developments in Computational Methods at the Nanoscale[END_REF][START_REF] Heinz | Simulations of inorganic-bioorganic interfaces to discover new materials: insights, comparisons to experiment, challenges, and opportunities[END_REF]. Apart from the atomistic structural data necessary to construct the molecular models, the most fundamental information required to perform such simulations is contained in the sets of parameters describing interatomic interactions in the modeled systems, often collectively called the force fields. The fact that in the literature there are currently several widely used and well tested classical force fields available for molecular modeling of clays and related inorganic materials (e.g. [START_REF] Heinz | Force field for mica-type silicates and dynamics of octadecylammonium chains grafted to montmorillonite[END_REF] and even a larger number of various force fields for molecular simulations of organic and bio-organic molecules (e.g. [START_REF] Guvench | Comparison of protein force fields for molecular dynamics simulations[END_REF], often leaves unanswered the important question of finding an optimal combination of the force field parameters for accurate molecular modeling of both the organic and inorganic parts of the composite systems. Indeed, estimating the actual predictive capabilities of such simulations is a highly non-trivial task. In order to quantitatively address this issue one needs a testing model which should be known to provide a stable interlayer complex, and should be supported by sufficient amount of reliable experimental data for comparisons with the simulated results. For these reasons, complexes of smectites with ethylene glycol (EG)-water mixtures were selected for the present study.

Intercalation of EG in hydrated divalent ion-smectites is known to provide a structure with relatively constant basal spacing (e.g. [START_REF] Mosser-Ruck | Effects of ethylene glycol saturation protocols on XRD patterns: A critical review and discussion[END_REF]. During this process, EG molecules penetrate into the interlayer spaces of the swelling clays, leading to the formation of a two-layer structure (d ~ 17 Å). For the purpose of X-ray diffraction data interpretation, a simplified model of this complex has been proposed by [START_REF] Reynolds | An X-ray study of an ethylene glycol-montmorillonite complex[END_REF].

Further studies have shown that the basal spacing is larger for clay minerals with lower layer charge, but it also depends on the localization of the charge, the type of the exchangeable cations, the particle size and the relative humidity [START_REF] Harward | Swelling properties of synthetic smectite in relation to lattice substitutions[END_REF][START_REF] Brindley | Ethylene glycol and glycerol complexes of smectite and vermiculites[END_REF][START_REF] Harward | Properties of vermiculite and smectites: Expansion and collapse[END_REF][START_REF] Środoń | Precise identification of illite/smectite interstratification by X-ray powder diffraction[END_REF][START_REF] Sato | Effects of layer charge, charge location, and energy change on expansion properties of dioctahedral smectites[END_REF]. If the charge is located in the tetrahedral sheet then the values of basal spacing are lower than those observed when the charge is located in the octahedral sheet [START_REF] Sato | Effects of layer charge, charge location, and energy change on expansion properties of dioctahedral smectites[END_REF]. Later studies have shown that smectites may form a one-layer EG complex instead of the two-layer complex at very low relative humidities [START_REF] Eberl | Sericite from the Silverton caldera, Colorado: Correlation among structure, composition, origin, and particle thickness[END_REF], or with K + as the exchange cation even at intermediate humidities [START_REF] Eberl | Potassium fixation in smectite by wetting and drying[END_REF]. The experimental observations were to some extent explained with the help of molecular simulations, e.g. the preference of formation of bilayer structure and the preferred number of water molecules in the structure for the Ca 2+ form [START_REF] Szczerba | Ethylene glycol intercalation in smectites. Molecular dynamics simulation studies[END_REF]. Several methods of clay glycolation with specific technicalities are possible and [START_REF] Mosser-Ruck | Effects of ethylene glycol saturation protocols on XRD patterns: A critical review and discussion[END_REF] have demonstrated certain differences and inconsistencies between them. The formation of monolayer or bilayer complexes was found to be dependent on the type of the interlayer cation and the concentration of EG in the close vicinity around smectite (the glycolation protocol). There is also a gradual loss of EG from the complex with smectite because of the equilibration between the complex and the surrounding environment.

Therefore, X-ray measurements on such samples should be performed relatively quickly, within a few hours after glycolation. The dynamics of EG intercalation into smectites was studied by [START_REF] Svensson | Intercalation of smectite with liquid ethylene glycol -Resolved in time and space by synchrotron X-ray diffraction[END_REF] using synchrotron X-ray diffraction techniques. This study showed that EG molecules replace H2O in the interlayer, however, still leaving significant amounts of water in the structure.

The present study was undertaken in order to investigate the effects of different combinations of force field parameters on the accuracy of two-EG layer structure representation. For this purpose inorganic (clay) and organic (EG) force field parameters need to be coupled in one simulation.

For clay minerals there exist several force field parameterizations that were used for molecular simulations of these structures over the past three decades. These parameters can be preliminarily divided into two groups: the ones that rely on the assumption of a rigid clay framework and the ones allowing partial orfull flexibility of the structure. Historically, the parametrizations of the first kind were developed earlier (e.g. [START_REF] Skipper | Computer simulation of interlayer water in 2:1 clays[END_REF][START_REF] Smith | Molecular computer simulations of the swelling properties and interlayer structure of cesium montmorillonite[END_REF]. They are, however, inherently limited and can lead to incorrect representation of the processes of adsorption, surface hydration, hydrogen bonding, diffusion rates etc. (e.g. [START_REF] Cygan | Molecular models and simulation of layered minerals[END_REF]. The flexible force field parameterizations are currently much wideer represented (e.g. [START_REF] Hill | Molecular mechanics potential for silica and zeolite catalysts based on ab initio calculations. 2. Aluminosilicates[END_REF][START_REF] Teppen | Molecular dynamics modeling of clay minerals. 1. Gibbsite, kaolinite, pyrophillite, and beidellite[END_REF][START_REF] Sato | Molecular simulation for flexibility of a single clay layer[END_REF][START_REF] Manevitch | Elastic properties of a single lamella of montmorillonite by molecular dynamics simulation[END_REF]Ruthlege, 2004, Cygan et al., 2004;[START_REF] Heinz | Force field for mica-type silicates and dynamics of octadecylammonium chains grafted to montmorillonite[END_REF]. It was shown, however, that only the parametrizations of [START_REF] Heinz | Force field for mica-type silicates and dynamics of octadecylammonium chains grafted to montmorillonite[END_REF] and ClayFF [START_REF] Cygan | Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field[END_REF] give cleavage energies that reasonably correspond to the experimental values [START_REF] Heinz | Force field for mica-type silicates and dynamics of octadecylammonium chains grafted to montmorillonite[END_REF]. Therefore, only these force fields were tested in this study. Additionally, [START_REF] Ferrage | Hydration Properties and Interlayer Organization of Water and Ions in Synthetic Na-Smectite with Tetrahedral Layer Charge. Part 2. Toward a Precise Coupling between Molecular Simulations and Diffraction Data[END_REF] have shown that a better description of the hydrated interlayer structure with the ClayFF force field can be achieved if the Lennard-Jones parameters of clay surface oxygens are increased by ~7%.

The interatomic interaction parameters for EG can be selected from a wide range of different organic force fields. Usually, these force fields are parameterized to optimize the description of certain specific sets of organic molecules. For example, CHARMM was constructed to model proteins (MacKerrel et al., 1998), AMBER -to model peptides, proteins and nucleic acids [START_REF] Cornell | A second generation force field for the simulation of proteins, nucleic acids and organic molecules[END_REF], OPLS-aa for simulations of liquid hydrocarbons [START_REF] Jorgensen | Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids[END_REF] and CVFF for amino acids, hydrocarbons, and many other organic molecules [START_REF] Dauber-Osguthorpe | Structure and energetics of ligand binding to proteins: E. coli dihydrofolate reductase-trimethoprim, a drug-receptor system[END_REF]. More recent parameterizations which belong to a certain earlier family but were later re-optimized for a larger set of organic molecules include CHARMM-related CGenFF (general force field for drug-like molecules; [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] and AMBER-related GAFF (general AMBER force field; [START_REF] Wang | Development and testing of a general amber force field[END_REF]. For the purpose of this study the parametrization of a given family that was best optimized for small organic molecules, such as EG, was taken into account.

Even a more diverse set of force field parameterizations for water molecules are available in the literature (e.g. [START_REF] Wallqvist | Molecular Models of Water: Derivation and Description[END_REF][START_REF] Guillot | A reappraisal of what we have learnt during three decades of computer simulations on water[END_REF][START_REF] Kalinichev | Molecular simulations of liquid and supercritical water: Thermodynamics, structure, and hydrogen bonding[END_REF]. In this study only the SPC (simple point charge)model of [START_REF] Berendsen | Interaction models for water in relation to protein hydration[END_REF] is considered because it is known to perform well in combination with ClayFF, gives relatively good solution structure for smectite interfaces and interlayers (e.g. [START_REF] Heinz | Force field for mica-type silicates and dynamics of octadecylammonium chains grafted to montmorillonite[END_REF][START_REF] Morrow | Structure, energetics, and dynamics of smectite clay interlayer hydration: Molecular dynamics and metadynamics investigation of Na-hectorite[END_REF][START_REF] Ngouana Wakou | Structural arrangements of isomorphic substitutions in smectites: Molecular simulation of the swelling properties, interlayer structure, and dynamics of hydrated Cs-montmorillonite revisited with new clay models[END_REF][START_REF] Greathouse | Molecular simulation of structure and diffusion at smectite-water interfaces: Using expanded clay interlayers as model nanopores[END_REF] and is also consistent with the CVFF and INTERFACE parametrizations. Thoretically, the application of different force field parametrizations for H2O molecules can affect the simulation results for hydrated clays. The work of [START_REF] Ferrage | Hydration Properties and Interlayer Organization of Water and Ions in Synthetic Na-Smectite with Tetrahedral Layer Charge. Part 2. Toward a Precise Coupling between Molecular Simulations and Diffraction Data[END_REF] has shown, however, that using the SPC model is sensible in terms of the agreement of MD simulation results with X-ray diffraction data.

Hybrid organo-clay materials is a rapidly growing area of molecular simulation studies that requires a reliable coupling of different force field parameters. Very often ClayFF is used together with CVFF (e.g. [START_REF] Kumar | Hydration, swelling, interlayer structure, and hydrogen bonding in organolayered double hydroxides: Insights from molecular dynamics simulation of citrate-intercalated hydrotalcite[END_REF][START_REF] Liu | Interlayer structure and dynamics of alkylammonium-intercalated smectites with and without water: A molecular dynamics study[END_REF][START_REF] Suter | Computer simulation study of the materials properties of intercalated and exfoliated poly (ethylene) glycol clay nanocomposites[END_REF][START_REF] Kalinichev | Molecular dynamics computer simulations of the effects of hydrogen bonding on the properties of layered double hydroxides intercalated with organic acids[END_REF]. There are also examples of successful use of ClayFF in combination with OPLS-aa (e.g. [START_REF] Schampera | Surface structure of organoclays as examined by X-ray photoelectron spectroscopy and molecular dynamics simulations[END_REF], CHARMM (e.g. [START_REF] Duque-Redondo | Molecular Forces Governing Shear and Tensile Failure in Clay-Dye Hybrid Materials[END_REF] and AMBER (e.g. [START_REF] Wang | Molecular adhesion at clay nanocomposite interfaces depends on counterion hydration-Molecular dynamics simulation of montmorillonite/xyloglucan[END_REF][START_REF] Swadling | Clay minerals mediate folding and regioselective interactions of RNA: a large-scale atomistic simulation study[END_REF]. INTERFACE force field [START_REF] Heinz | Thermodynamically consistent force fields for the assembly of inorganic, organic, and biological nanostructures: The INTERFACE force field[END_REF] is also recently developed to be used in the studies of interactions of organic molecules with minerals, including smectite clays, together with CVFF, PCFF and COMPASS organic parameterizations. There are also other examples of combinations of older clay mineral force field parameters with different organic force fields: e.g. [START_REF] Tambach | Hysteresis in clay swelling induced by hydrogen bonding: accurate prediction of swelling states[END_REF] used force field of [START_REF] Skipper | Monte Carlo simulation of interlayer molecular structure in swelling clay minerals. I: Methodology[END_REF] together with OPLS-ua [START_REF] Jorgensen | Monte Carlo simulations of the hydration of ammonium and carboxylate ions[END_REF]. Not all of these earlier force field combinations were extensively tested before they we used in this area of very active current research.

At the same time, the structures obtained from molecular simulations using different force field parameters can be validated by their comparison with available X-ray diffraction data. This methodology has been widely used for hydrated clay minerals (e.g. [START_REF] Ferrage | Hydration Properties and Interlayer Organization of Water and Ions in Synthetic Na-Smectite with Tetrahedral Layer Charge. Part 2. Toward a Precise Coupling between Molecular Simulations and Diffraction Data[END_REF]. Even though only the atomic density distribution of the interlayer species along a direction perpendicular to the layering plane can be considered for turbostratic smectites, it is still a very powerful methodology. It was used in the present study for the purpose of finding the optimal set of force field parameters for organic molecules interacting with smectites.

METHODOLOGY

Smectite samples used in the study

Three smectite samples from the collection of the Source Clays Repository of The Clay Minerals Society were studied. They have crystal structures with notably different layer charge and charge localization: low charge montmorillonite (SWy-1), high charge montmorillonite (SAz-1), and high charge beidellite (SbCa-1). Some of their characteristics are summarized in Table 1.

=== Table 1 ===

The particle size fractions were separated to avoid contamination with phases other than smectite. All samples were studied in Ca 2+ form that was prepared via dialysis.

The samples were sedimented from aqueous solution on a glass slide. Then all slides were glycolated with the following procedure: each air-dried sample from RH around 50% was put into a closed plastic box (15×15×3 cm) that was placed in a vessel with EG solution poured at the bottom. The sample within the box was heated at 60 o C for 10 h in an oven.

XRD patterns registration

To minimize EG loss from the prepared EG-smectite complexes the XRD patterns for all samples were measured directly after glycolation and in EG-saturated atmosphere. The Xray diffraction datawere collected in the range from 2 o to 60 o 2Θ using a Thermo ARL XRD system, CuKα radiation, and a Peltier-cooled solid-state detector. The tube current and voltage were 45 mA and 35 kV, respectively. The following slit sizes from tube to detector were used: 0.9 mm (0.645 o ), 1.3 o Soller, 1.05 mm, sample, 1.0 mm, 1.3 o Soller, 0.3 mm. The step size was 0.05 o and the counting time was 20 s per step.

Molecular dynamics simulations

The simulated structural models of smectites were based on the pyrophyllite crystallographic data [START_REF] Lee | Single crystal X-ray refinement of pyrophyllite-1Tc[END_REF], with several isomorphic substitutions introduced at particular atomic sites to mimic the experimentally studied smectite samples. The total energy of a molecular model is usually described by a sum of Coulombic (electrostatic) interactions, short-range non-electrostatic interactions (sometimes referred to as the Van der Waals terms), and bonded (intramolecular) interactions:

Bonded VDW Coul Total E E E E + + = (1)
The bonded terms are especially important for organic molecules and typically include the bond stretching, angle bending energy terms, various torsional terms, etc. The electrostatic energy is represented by Coulomb's law:

∑ ≠ = j i ij j i Coul r q q e E o 2 4πε , ( 2
)
where rij is the separation distance between the charged atoms i and j, e is the charge of the electron, εo is the dielectric permittivity of vacuum (8.85419 × 10 -12 F/m), and the partial charges qi and qj are usually derived from quantum mechanics calculations and assigned by the specific force field model. The Van der Waals energy term is usually represented by the conventional Lennard-Jones (12-6) function, and includes the short-range repulsion associated with the increase in energy as two atoms closely approach each other and the attractive dispersion energy:

∑ ≠                 -         = j i ij ij o ij ij o ij o VDW r R r R D E 6 , 12 , , 2 .
(3) Do,ij and Ro,ij are empirical parameters specific to a particular force field model. The interaction parameters between the unlike atoms are usually calculated according to the socalled Lorentz-Berthelot mixing rules (e.g. [START_REF] Allen | Computer Simulation of Liquids[END_REF]: arithmetic mean for the distance parameter, Ro, and the geometric mean for the energy parameter, Do:

) R R ( R j , o i , o ij , o + = 2 1 (4) 10 j , o i , o ij , o D D D = (5)
The Coulombic and VDW interactions are excluded for proximate intramolecular (bonded) interactions (i.e., 1-2 and 1-3 atom position exclusions) when large organic molecules are modeled.

In our case, the interatomic interactions involving smectite structures were described with the ClayFF [START_REF] Cygan | Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field[END_REF] and the INTERFACE [START_REF] Heinz | Force field for mica-type silicates and dynamics of octadecylammonium chains grafted to montmorillonite[END_REF] force fields.

GAFF, OPLS-aa, CGenFF, and CVFF organic force fields were used for EG, while the SPC model was assumed for water molecules [START_REF] Berendsen | Interaction models for water in relation to protein hydration[END_REF]. The INTERFACE force field was only coupled with the CVFF organic force field because of the consistency of its 1-4-scaling factors (value of 1.0). CGenFF assumes no scaling, while OPLS-aa uses a scaling factor of 0.5, and GAFF uses 0.5 for the non-bonded non-electrostatic part of the interaction and 0.8333 for the electrostatic interactions. ClayFF has no restrictions on the organic force field that it can be coupled with, because no 1-4 pairs are present in this parametrization.

Additionally, in the original clay mineral force fields, the Lennard-Jones parameters for basal oxygens, Ob, were modified as suggested by Ferrage et al. (2011) (ClayFFmod and INTERFACEmod). The LJ parameters of tetrahedral Si and Al atoms were, respectively, also modified in these models to maintain the same Si-Ob and Al-Ob distances as in the original force fields (see [START_REF] Szczerba | Structure and dynamics of water-smectite interfaces: Hydrogen bonding and the origin of the sharp O-Dw/O-Hw infrared band from molecular simulations[END_REF] and supplementary materials for a more detailed discussion). Ewald summation was applied to calculate the long range corrections to the Coulombic interactions (e.g., [START_REF] Allen | Computer Simulation of Liquids[END_REF] and the cut-off distance was set to 8.5 Å.

A set of different EG compositions around those of Reynolds (1965) -1.7 EG per half unit cell (phuc) -was covered by the simulations to obtain the basal spacing closest to the experimentally observable. A constant number of interlayer water molecules -0.8 H2O phuc -was assumed. Additionally, for the set of parameters that gave the best results in previous calculations, the EG content was set to vary from 1.4 to 2.0 with a step of 0.2 phuc, while the water content also varied between 0.0 and 1.2 H2O phuc with a step of 0.3 phuc.

For all models, structural optimizations (total energy minimizations) were performed first, followed by NPT-ensemble MD simulations at 1 bar under evolving temperature using a Langevin dynamics algorithm to control the temperature and a Langevin piston to control the pressure. The time step to integrate the equations of atomic motion was set to 1 fs and the dynamic trajectories of all atoms and system properties were recorded every 1 ps. For the first 0.5 ns the temperature was set to 398 K. Then the temperature dropped to 298 K and the simulation continued for 1 ns. From the last 0.5 ns the equilibrium system properties were recorded for further analysis. To exclude any undesirable displacement of the center of mass of the simulated model, in all simulations one atom of the octahedral sheet was fixed at its initial position, but its interactions with all neighboring atoms were still fully accounted for.

All MD simulations were performed using the LAMMPS computer program [START_REF] Plimpton | Fast parallel algorithms for short-range molecular dynamics[END_REF] http://lammps.sandia.gov).

Parametrization of EG molecules with organic force fields

CGenFF parameters were automatically assigned with the help of the on-line tool (cgenff.paramchem.org) using the NAMD convention, which was modified to make it consistent with the LAMMPS format. GAFF parametrization was performed using the moltemplate program (moltemplate.org) with charges automatically assigned using the program TPACM4 available at www.scfbio-iitd.res.in/software/drugdesign/charge.jsp [START_REF] Mukherjee | A Fast empirical GAFF compatible partial atomic charge assignment scheme for modeling interactions of small molecules with biomolecular targets (TPACM4)[END_REF]. CVFF parameters were automatically assigned within this force field's implementation in the LAMMPS program. Partial atomic charges were calculated based on the bond increments from the file containing CVFF parameters. OPLS-aa parameters were also automatically assigned using the moltemplate program. Tcl scripts executed in the VMD molecular visualization program [START_REF] Humphrey | VMD -Visual Molecular Dynamics[END_REF] were used to generate all input files for LAMMPS simulations.

Calculation of the simulated XRD patterns

For every simulated structure, the atomic density distributions of the interlayer species as well as electron density profiles along the direction perpendicular to the layering were calculated. Then all the distributions were symmetrized with respect to the interlayer center and used as an input for calculations of the X-ray diffraction patterns using the Sybilla code (Chevron proprietary). The values of Tmean, d-spacing, σ* and ∆-d-spacing were optimized automatically. The range of 2θ < 4.5 o was excluded from the optimization because in this range the effect of super-crystallites becomes important, which is not taken into account by the code. The amount of Fe in the octahedral sheet was taken directly from the Table 1 and was not optimized.

Tmean is the average thickness calculated for a lognormal distribution of crystallite thicknesses. A variation of this value affects mainly the broadening of the 00l peaks.The dspacing represents the layer-to-layer distance of smectite structures and principally affects only the positions of 00l maxima. The parameter of σ* is a standard deviation of the Gaussian orientation function of crystallites [START_REF] Reynolds | The Lorentz-polarization factor and preferred orientation in oriented clay aggregates[END_REF]. The ∆-d-spacing describes fluctuations of the layer-to-layer distance. The variations of σ* and ∆-d-spacing parameters affect the relative XRD intensities by a factor that is quite complex to express analytically, but is known to be monotonic function of 2θ. Generally, it can be concluded that the largest modification of the relative intensities of 00l reflections is influenced by the structure of interlayer species, which depends on the force-fields used. This is also because all the parameters were optimized to reproduce very similar values for a certain studied smectite.

RESULTS AND DISCUSSION

Comparison of organic force fields for EG

In order to compare the performance of organic force fields in reproducing bulk liquid EG properties, the density of a box consisting of 230 EG molecules as well as the angular distributions of the torsional O-C-C-O potential energy terms were calculated. Table 2 compares the calculated liquid densities with the experimental data. The closest agreement is achieved by the GAFF parameterization, while all other force fields result in densities noticeably lower (OPLS-aa and CGenFF) or higher (CVFF) than the experimental value.

=== Table 2 ===

The torsional O-C-C-O angle potential energy distributions were compared to a reference obtained from DFT calculations performed under the B3LYP/DGDZVP level of theory using Gaussian Inc. software package [START_REF] Frisch | Gaussian-94[END_REF]. In these calculations, all distances and angles were optimized for a certain O-C-C-O angle. The results (Figure 1) demonstrate that the CVFF and OPLS-aa force field parameterizations are the closest to the reference quantum chemical calculation. In the case of GAFF parameterization, the gauche conformation is predicted to be too stable relatively to the trans one. On the other hand, the stability of the trans conformation is overestimated in the CGenFF parameterization. All these differences are important for the proper modeling of EG conformations in smectite interlayers.

=== Figure 1 ===

Comparison of intercalate structures with 1.7 EG and 0.8 H2O phuc

The observed differences between the organic force fields, described in the previous section, should have an effect on the simulated interlayer structure and resulting X-ray diffraction patterns. To study this, the intercalate structure with EG and water content corresponding to that of [START_REF] Reynolds | An X-ray study of an ethylene glycol-montmorillonite complex[END_REF] (1.7 EG and 0.8 H2O phuc) was calculated first. An example for SWy-1 is shown in Figure 2 (the results for two other smectites are provided in the supplementary materials).

=== Figure 2 ===

The structures obtained using the INTERFACE and INTERFACEmod force fields are substantially different than those resulting from the application of the ClayFF and ClayFFmod -i.e., oxygen atoms of EG molecules have much higher tendency to be located closer to the surface than their carbon atoms. In the case of the latter force fields, the positions of the closest to the surface carbon and oxygen atoms do not significantly differ. A similar effect is also observed for water molecules: they are located closer to the clay surface for the INTERFACE-than for the ClayFF-based models.

Different organic force fields have also substantial effect on the resulting averaged structures of the intercalate. The effect of different EG liquid density is reflected in the differences in the basal spacings -the lower the density, the higher the basal spacing. There is no obvious relation between the observed relative stabilities of the trans-gauche conformers (Figure 1) and the obtained interlayer structures. The position of carbon atoms is relatively similar for all the organic force fields tested, while there are substantial differences in the positions of the oxygen atoms of EG. These differences, however, do not affect substantially the calculated electron density profiles (Figure 3 and supplementary materials).

The modification of the Lennard-Jones parameters for the basal oxygens of smectites leads to a very slight modification of the basal spacing. However, it affects the interlayer structure in such a way that organic and H2O molecules are located at somewhat larger distances from the clay surface than in the original ClayFF models (Figure 2). The same tendency was observed by [START_REF] Ferrage | Hydration Properties and Interlayer Organization of Water and Ions in Synthetic Na-Smectite with Tetrahedral Layer Charge. Part 2. Toward a Precise Coupling between Molecular Simulations and Diffraction Data[END_REF]. These changes, in turn, affect the calculated electron density profiles (Figure 3).

=== Figure 3 ===

Comparison of intercalate structures with 0.8 H2O phuc and basal spacing close to experimental

Based on the structure proposed by [START_REF] Reynolds | An X-ray study of an ethylene glycol-montmorillonite complex[END_REF] and on the previous MD simulation studies of EG-smectite structure [START_REF] Szczerba | Ethylene glycol intercalation in smectites. Molecular dynamics simulation studies[END_REF], the content of water is estimated to be around 0.8 -1.0 molecules phuc. This value was therefore assumed as a constant and chosen to be equal 0.8 phuc in the present simulations. The EG content was adjusted to achieve the basal spacing close to the experimental value. This is an obvious approximation due to the fact that the content of water in the structure will be related to the actual relative humidity and the time between the sample glycolation and XRD measurement and the conditions of the XRD profile registration. The content of water and EG is also dependent on the glycolation procedure [START_REF] Mosser-Ruck | Effects of ethylene glycol saturation protocols on XRD patterns: A critical review and discussion[END_REF]. Without this approximation, however, the number or necessary calculations would increase dramatically, but without affecting the main conclusions too much.

The results of the calculations for SWy-1 are presented in Figure 4 (the results for two other smectites are presented in the supplementary materials). The effect of the variation in the EG content on the structures is not very strong (in comparison to Figure 2), and the main factors affecting the distribution of atoms in the interlayer space are the organic force field as well as the clay mineral force field.

=== Figure 4 ===

Comparison of the calculated X-ray diffractograms

In order to evaluate and further quantify the differences due to the application of different organic force fields in the molecular modeling studies of organic-smectite interactions, all the above calculated distributions were incorporated into the Sybilla code to calculate X-ray diffractograms for the simulated structures. Thus calculated diffractograms with optimized values of Tmean, d-spacing, σ*, and ∆-d-spacing were then compared with experimental ones (Figures 567). The results show quite significant discrepancies in the obtained XRD patterns for different layer charges as well as for different combinations of force fields. Because XRD is only sensitive to the distribution of electrons (the simulated electron density profiles are shown in Figure 3), different electron density profiles may results in different calculated XRD patterns, especially in the distribution of relative 00l intensities.

For low-charge montmorillonite (SWy-1; charge 0.27 phuc) the calculated results agree well with experiment for all clay mineral force fields, and only relatively small differences between the organic force fields used are noticeable. The modification of the Lennard-Jones parameters for the basal oxygens leads to some increase of intensities of the 002 and 003 peaks and to a small decrease of intensities of the 004, 005 and 006 peaks.

For high-charge montmorillonite (SAz-1; charge 0.54 phuc) and beidellite (SbCa-1; charge 0.50 phuc), the modification of the Lennard-Jones basal oxygens parameters of the ClayFF-based models substantially improves the agreement between the theoretical and experimental diffractograms for all the organic force fields tested. A substantial correction of the 002 and 003 peak intensities towards those of experimental values is observed. A similar modification applied to the INTERFACE-based models leads to only relatively small improvement. This indicates that, although this force field provides a good structure for lowcharge montmorillonites, it should be probably modified to describe more accurately the interlayer structure of high-charge smectites.

Among the organic force fields, the least accurate results are obtained with CGenFF. This is probably related to its overestimation of the trans EG molecular conformation and its poor agreement of the EG liquid density with experiment (Figure 1 and Table 2). Combined with unmodified ClayFF, GAFF gives the best results. The resulting X-ray diffractogram is, however, poor for high charge-beidellite in any combination of organic force field with unmodified clay mineral set of parameters. The two other sets (OPLS-aa and CVFF) give the best results for ClayFFmod. The observed differences can be related to the approximations of the organic, clay mineral and water force fields and possible uncertainties in the number of EG and water molecules in the interlayer estimated from experimental data and assumed in the models.

=== Figure 5 === === Figure 6 === === Figure 7 ===

Comparison of the structures with different water and EG contents

Based on the results presented in the previous section, one of the two best sets of force field combinations (ClayFFmod + CVFF) was used to study the interlayer structures further and to check if the variation of the assumed EG and water content can improve the calculated X-ray diffractograms. A range of EG compositions was systematically probed between 1.4 and 2.0 phuc with a step of 0.2 phuc, with water contents varying between 0.0 and 1.2 phuc with a step of 0.3 phuc. The variation of d-spacing depending on the EG and water content (Figure 8) shows a plateau at around 16.5 -17.0 Å with some further variation depending on the layer charge. Based on these plots, several structures were selected (black circles in Figure 8) to calculate average atomic distributions and then to calculate the corresponding X-ray diffraction patterns for two-layer intercalate structures.

The best of the calculated X-ray diffractograms shows only a small improvement compared to the results shown in figures 5-7 (Figure 9; all results are presented in the supplementary materials). The discrepancy between the simulated and experimental diffractograms can, therefore, be primarily attributed to the inaccuracy of the force fields used. The most important factor is clearly the selection of the clay mineral force field -its parameters may require additional improvement.

=== Figure 8 === === Figure 9 ===

CONCLUSIONS

The main factor affecting the accuracy of the calculated X-ray diffractograms of EGwater intercalate in smectites is the force field selected to model clay mineral substrate. The selection of organic force field parameters, although also important, has only a secondary effect on the results obtained. The best set of parameters for smectite was found to be ClayFF with modified Lennard-Jones parameters of the basal surface oxygens (ClayFFmod). All tested organic force fields perform relatively well in combination with ClayFFmod and small differences between them depend on the value and location of the smectite charge. Generally, among the organic force fields tested, GAFF, OPLS-aa, and CVFF perform relatively well, while the application of CgenFF leads to the least accurate results.

The INTERFACE force field produces relatively good results for low-charge montmorillonite, but it performs less well for high-charge smectites. Unlike for ClayFF, modifications of the Lennard-Jones parameters for INTERFACE do not help much to improve the resulting simulated X-ray diffractograms. The origin of the discrepancy between the simulated and experimental diffractograms is apparently in the overestimation of the interlayer atomic populations close to the clay surface.

Further assessment of molecular models could be performed using neutron diffraction. This would, in particular, help in discriminating between models giving relatively similar electron density profiles but contrasting distributions of interlayer hydrogen atoms. Using this approach, positions of water molecules in the interlayer space of smectite complexes with organic molecules can be much more precisely determined.

  All three structures were built by substituting a relevant number of Al atoms with Mg and Si with Al in the octahedral and tetrahedral sheets of the clay structure, respectively. The Mg/Al ordering in the octahedral sheets was introduced following the work of Ortega-Castro et al. (2010), i.e. maximizing the distance between Mg atoms. The Al/Si ordering in the tetrahedral sheet was random but obeying the Löwenstein rule, i.e. exluding Al-O-Al linkages. The simulation supercell was 8×4×2 unit cells in the a, b, and c crystallographic directions, respectively (~41.6 Å × 36.1 Å × Z Å; the value of Z varied depending on the amount of EG and H2O in the interlayer space and the force field used).

Figure 1 .

 1 Figure 1. Potential energy scan of the torsional O-C-C-O angle of the EG molecule as reproduced by different organic force fields in comparison with quantum chemical calculations at the B3LYP/DGDZVP level of theory.

Figure 2 .

 2 Figure2. Z-density profile distributions for carbon and oxygen of EG, oxygen and hydrogen of water, and Ca 2+ ions for SWy-1 for 1.7 EG and 0.8 H2O phuc. The resulting basal spacings are also shown.

Figure 3 .

 3 Figure 3. Electron density profiles obtained from MD simulations (black lines) compared to the profiles suggested by Reynolds (1965) (grey lines) for SWy-1 for 1.7 EG and 0.8 H2O phuc.

Figure 4 .

 4 Figure 4. Z-density profile distributions for carbon and oxygen atoms of EG, oxygen and hydrogen atoms of water and Ca 2+ ions for SWy-1 with 0.8 H2O phuc and variable EG content for which the basal spacing is close to 16.92 Å. The resulting basal spacings are also shown.

Figure 5 .

 5 Figure 5. Comparison of SWy-1 XRD patterns for various force field combinations (red -experimental, blackcalculated).

Figure 6 .

 6 Figure 6. Comparison of SAz-1 XRD patterns for various force field combinations (red -experimental, blackcalculated).

Figure 7 .

 7 Figure 7. Comparison of SbCa-1 XRD patterns for various force field combinations (red -experimental, black -calculated).

Figure 8 .

 8 Figure 8. d-spacing dependence on the EG and water content for: a) SWy-1, b) SAz-1, c) SbCa-1 smectites.Small circles mark the structures which were used to calculate X-ray diffraction patterns. Thick black lines correspond to experimental basal spacings.

Figure 9 .

 9 Figure 9. Comparison of XRD patterns for different smectites selected from the compositions in figure 8, that provide the best agreement between experimental and theoretical diffractograms.

Fig. 1

 1 Fig.1

Fig. 2

 2 Fig.2

Fig. 3

 3 Fig.3

Fig. 4

 4 Fig.4

Fig. 5

 5 Fig.5

Fig. 6

 6 Fig.6

Fig. 7

 7 Fig.7

Fig. 8

 8 Fig.8

Fig. 9

 9 Fig.9

Table 1 .

 1 Smectites used in the study.

					Particle size
	Smectite	Tetrahedral charge	Octahedral charge	Fe in octahedral sheet
					fraction
	SWy-1	0.0	0.28	0.20	<0.1 μm
	SAz-1	0.0	0.56	0.26	<2.0 μm
	SbCa-1	0.50	0.0	0.09	<1.0 μm

Table 2 .

 2 Comparison of calculated and experimental liquid EG density values

	Organic force field	Density at 20 o C (g/cm 3 )
	OPLS-aa	1.060
	GAFF	1.123
	CGenFF	1.054
	CVFF	1.188
	Experimental	1.115 a

a The experimental value of EG density at 20 o C is taken from Dow Chemicals (dow.com/ethyleneglycol/about/properties.htm).

ACKNOWLEDGMENTS

This project was made possible with financial support from National Science Centre (grant 2012/05/B/ST10/01948), which provided IGS PAS with high-performance computing server.

Computer time allocations made available within the Distributed European Computing

Initiative (project DEC07_NUWCLAY and DEC11 COMPCLAY by the PRACE-2IP receiving funding from the European Community's FP7/2007-2013 under grant agreement RI-283493) and within PLGRID infrastructure are also gratefully acknowledged.

AGK was also supported by the industrial chair "Storage and Disposal of Radioactive Waste" at the Ecole des Mines de Nantes, funded by ANDRA, Areva, and EDF.