Pier Paolo Poier
email: pier.poier@sorbonne-universite.fr

Théo Jaffrelot Inizan

Olivier Adjoua

Louis Lagardère

Jean-Philip Piquemal
email: jean-philip.piquemal@sorbonne-universite.fr

Accurate Deep Learning-aided Density-free Strategy for Many-Body Dispersion-corrected Density Functional Theory

Using a Deep Neuronal Network model (DNN) trained on the large ANI-1 data set of small organic molecules, we propose a transferable density-free many-body dispersion model (DNN-MBD). The DNN strategy bypasses the explicit Hirshfeld partitioning of the Kohn-Sham electron density required by MBD models to obtain the atomin-molecules volumes used by the Tkatchenko-Scheffler polarizability rescaling. The resulting DNN-MBD model is trained with minimal basis iterative Stockholder atomic volumes and, coupled to Density Functional Theory (DFT), exhibits comparable (if not greater) accuracy to other approaches based on different partitioning schemes. Implemented in the Tinker-HP package, the DNN-MBD model decreases the overall computational cost compared to MBD models where the explicit density partitioning is performed. Its coupling with the recently introduced Stochastic formulation of the

Since its original formulation in 1965, Kohn-Sham Density Functional Theory 1 (KS-DFT) has become the most popular family of electronic structure methods. KS-DFT represents in fact the cheapest way for introducing electronic correlation as its computational cost is similar to that of the Hartree-Fock method. KS-DFT is based on the idea of evaluating the kinetic energy from a Slater determinant thus assuming the electrons to be non-interacting. This apparently crude assumption actually leads to big improvements in describing chemical bonding compared to, for example, the use of the Thomas-Fermi kinetic energy formulation. The difference between the Slater determinant kinetic energy representation and the true one, together with the difference between the true total electronic interaction and the exchange energies represents, in KS-DFT, the key contribution to the exchange-correlation functional which remains, however, unknown.

In practice, the plethora of existing KS-DFT variants differentiate themselves in the way the exchange-correlation functional is approximated. Typically it is assumed to be a functional of the local electron density and eventually of its gradient and Laplacian. As a consequence, only local contributions to electronic correlation are included and this explains the general inadequacy of DFT methods to describe dispersion interactions which, on the other hand, have roots in long-range electronic correlation.

To retain the pleasant computational performances of KS-DFT methods, several dispersion corrections have been proposed. [START_REF] Grimme | Dispersion-corrected mean-field electronic structure methods[END_REF] Among these, the popular and successful approach of Grimme includes dispersion via empirical pairwise C 6 terms. [START_REF] Grimme | Accurate description of van der Waals complexes by density functional theory including empirical corrections[END_REF][START_REF] Grimme | Semiempirical GGA-type density functional constructed with a long-range dispersion correction[END_REF][START_REF] Grimme | A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu[END_REF] This is particularly appealing in virtue of its nearly zero additional computational cost.

A further approach is to replace the empirical pairwise terms with ones obtained from quantities coupled to the molecular electron density. For example, in Becke and Johnson's model, pairwise C 6 coefficients are written in terms of atomic polarizabilities and the averaged exchange-hole dipoles corresponding to each of the two atoms in the pair. [START_REF] Johnson | A post-Hartree-Fock model of intermolecular interactions: inclusion of higher-order corrections[END_REF][START_REF] Becke | Exchange-hole dipole moment and the dispersion interaction revisited[END_REF] In the alternative approach proposed by Tkatchenko and Scheffler (TS), [START_REF] Tkatchenko | Accurate molecular van der Waals interactions from ground-state electron density and free-atom reference data[END_REF] pairwise C 6 coefficients are instead expressed in terms of accurate free atom reference data as well as atoms-in-molecule (AIM) polarizabilities obtained from the rescaling of the corresponding AIM volumes computed via the Hirshfeld partitioning of the molecular electron density. [START_REF] Hirshfeld | Bonded-atom fragments for describing molecular charge densities[END_REF] One limitation of the above mentioned pairwise approaches is the impossibility of capturing non-additive many-body dispersion (MBD) effects, which inclusion has recently been shown important in modelling extended systems, supramolecular complexes and proteins in solutions, among others. [START_REF] Reilly | Seamless and accurate modeling of organic molecular materials[END_REF][START_REF] Ambrosetti | Hard numbers for large molecules: toward exact energetics for supramolecular systems[END_REF][START_REF] Ambrosetti | Wavelike charge density fluctuations and van der Waals interactions at the nanoscale[END_REF][START_REF] Stöhr | Quantum mechanics of proteins in explicit water: the role of plasmon-like solute-solvent interactions[END_REF] The non-additive long-range character of dispersion interactions has been modeled via a set of coupled fluctuating dipoles [START_REF] Langbein | Microscopic calculation of macroscopic dispersion energy[END_REF][START_REF] Donchev | Many-body effects of dispersion interaction[END_REF] (CFD) or alternatively by quantum Drude oscillators. [START_REF] Sommerfeld | Quantum Drude Oscillator Model for describing the interaction of rxcess electrons with water clusters: an application to (H2O)13[END_REF][START_REF] Jones | Quantum drude oscillators for accurate many-body intermolecular forces[END_REF][START_REF] Jones | Quantum Drude oscillator model of atoms and molecules: many-body polarization and dispersion interactions for atomistic simulation[END_REF][START_REF] Odbadrakh | Dispersion dipoles for coupled Drude oscillators[END_REF] In recent years, Tkatchenko, DiStasio Jr., Ambrosetti et al. have proposed a range-separated many-body dispersion model based on the CFD where the self-consistent screening of a set of atomic polarizabilities is performed (MBD@rsSCS). [START_REF] Tkatchenko | Accurate and efficient method for many-body van der Waals interactions[END_REF][START_REF] Ambrosetti | Long-range correlation energy calculated from coupled atomic response functions[END_REF] The MBD@rsSCS model is appealing not only for introducing non-additive many-body dispersion effects but also since it relies, de facto, on a single range-separation parameter which is tuned according to the choice of the exchange-correlation functional employed.

The MBD@rsSCS keeps in fact the spirit of the TS approach where AIM polarizabilities and van der Waals radii are obtained via the Hirshfeld partitioning of the density.

The Hirshfeld method leads to AIM densities which minimize the Kullback-Lieber divergence corresponding to the information loss upon molecule formation where this solid mathematical condition is used as a basis for the development of new information-theoretic partitioning methods. [START_REF] Heidar-Zadeh | Information-theoretic approaches to atoms-in-molecules: Hirshfeld family of partitioning schemes[END_REF] As discussed in references [START_REF] Ayers | Atoms in molecules, an axiomatic approach. I. Maximum transferability[END_REF][START_REF] Ayers | Variational principles for describing chemical reactions: condensed reactivity indices[END_REF] Hirshfeld partitioning makes its resulting AIM densities as close as possible to the ones of the isolated atoms, consequently AIM's properties turn out to be as similar as possible to those of the free atoms. This is particularly evident in the magnitude of Hirshfeld atomic charges, being too small in magnitude for reproducing the molecular electrostatic potential (ESP) or in modeling AIM polarizabilities in ionic and covalent crystals where the Hirshfeld partitioning leads to unrealistically large polarizabilities of cations which can even be found larger than those of the anions. [START_REF] Bučko | Extending the applicability of the Tkatchenko-Scheffler dispersion correction via iterative Hirshfeld partitioning[END_REF] The above mentioned shortcomings were ameliorated by the Iterative Hirshfeld (HI) scheme [START_REF] Bultinck | Critical analysis and extension of the Hirshfeld atoms in molecules[END_REF] where the reference atomic density employed in the partitioning is constructed as a linear combination of the two densities relative to the atomic oxidation states closest to the fractional number of electrons assigned by the partitioning at a given iteration.

The ESP computed from HI atomic charges have proven to agree remarkably well with ab initio computed reference. [START_REF] Van Damme | Electrostatic potentials from self-consistent Hirshfeld atomic charges[END_REF] In addition, the use of HI derived AIM polarizabilities leads to more realistic dispersion coefficients [START_REF] Bučko | Extending the applicability of the Tkatchenko-Scheffler dispersion correction via iterative Hirshfeld partitioning[END_REF] especially in ionic systems and adsorption phenomena on surfaces of ionic solids where the HI scheme used within the TS dispersion model improves interaction energies. 28 HI partitioning has also been employed in the MBD@rsSCS model replacing the original Hirshfeld scheme 29 and its use in in connection to the fractionally ionic AIM polarizabilities leads, in the just mentioned challenging systems, to reduced errors. [START_REF] Gould | A fractionally ionic approach to polarizability and van der Waals many-body dispersion calculations[END_REF] Despite the improvements carried by the HI partitioning, the scheme remains affected by a shortcoming arising from the density interpolation for negatively charged atoms as this procedure is, for some species, ill-defined. This arises from the fact that free anions such as N -and O 2 -(or in general any doubly negative ion) are not bound and their reference electron densities, computed at a complete basis set (CBS), result in a detached electron.

The iterative Stockholder atom (ISA) scheme, on the other hand, is not affected by this problem as the partitioning does not require reference atomic densities computed from isolated atoms at different ionic states as they are rather obtained from a spherical averaging of the molecular density using nuclei as expansion points. [START_REF] Lillestolen | Atomic charge densities generated using an iterative stockholder procedure[END_REF][START_REF] Misquitta | Distributed Multipoles from a Robust Basis-Space Implementation of the Iterated Stockholder Atoms Procedure[END_REF] The minimal basis iterative Stockholder atom (MBISA), a variant of the ISA method, have proven successful in the atomic polarizability rescaling approach employed by the TS scheme as well as in reproducing ab initio ESP from atomic point charges [START_REF] Verstraelen | Minimal basis iterative stockholder: atoms in molecules for force-field development[END_REF] and for this reasons its use in connection to the MBD@rsSCS model is particularly appealing. AIM properties are local quantities which depend on the near chemical environment and thus carry a certain degree of transferability. In particular, the TS polarizability rescaling scheme (employed in the MBD@rsSCS) makes use of AIM volumes which are well suited to be computed via deep neural network (DNN) where the environment vector associated to an atom's surrounding is defined within a local cutoff. The potential of deep learning in capturing local atomic properties has been proved by Isayev and co-workers whose multi-output DNN model successfully predicts AIM properties ranging from multipoles to volumes. [START_REF] Zubatyuk | Accurate and transferable multitask prediction of chemical properties with an atoms-in-molecules neural network[END_REF] In this Letter we present a hybrid DNN-aided MBD@rsSCS model (DNN-MBD) where the AIM volumes ratio employed in the TS polarizability rescaling are generated by a deep neural network trained on the ANI-1 data set (approximately 4.6 million structures) containing MBISA AIM volumes. [START_REF] Smith | The ANI-1ccx and ANI-1x data sets, coupled-cluster and density functional theory properties for molecules[END_REF] For the common S66x8 benchmark set, [START_REF] Řezáč | A Well-balanced database of benchmark interaction energies relevant to biomolecular structures[END_REF] the DNN-MBD model coupled to the common PBE/PBE0 density functionals, exhibits excellent interaction energies while completely bypassing the electron density partitioning with a consequent computational cost reduction. This electron density-free DNN-MBD approach employed in connection to our recently proposed linear scaling stochastic MBD@rsSCS formulation, [START_REF] Poier | N) stochastic evaluation of many-body van der Waals energies in large complex systems[END_REF] allows for modelling non-additive long-range dispersion interactions of up-to-millions atom systems at a very low computational cost without compromising the accuracy.

We note that kernel-ridge regression approaches to model AIM polarizabilities have been proposed in modelling dispersion interactions. [START_REF] Bereau | Non-covalent interactions across organic and biological subsets of chemical space: physics-based potentials parametrized from machine learning[END_REF][START_REF] Muhli | Machine learning force fields based on local parametrization of dispersion interactions: Application to the phase diagram of C 60[END_REF] This approach, however, is characterized by a O(N 2 ) and O(N 3 ) scaling of the required memory and computational cost involved in the model's training respectively, N being the size of the data set. These non-linear scaling prevents the applicability of kernel-ridge approaches on very large and diverse data sets, necessary for the generation of general-purpose MBD models. Additionally, the poor scaling with the number of processes limits its use on large systems. Here instead we generalize the approach to model MBD interactions to a much broader class of systems thanks to the employed model's flexibility and broad data set, without affecting the model's accuracy and linear scalability.

We will, in the following, proceed by briefly recalling the key concepts of the standard MBD@rsSCS model before introducing the DNN-MBD hybrid model and its performances.

As a starting point in this discussion, we examine the TS polarizability rescaling in Eq.( 1), where α i and V i represent the TS static polarizability and AIM volume respectively of the i-th atom while the zero superscript denotes free atom reference quantities.

α i = V i V 0 i α 0 i (1)
The AIM volume V i is obtained by solving the integral in Eq.( 2) where ρ(r) is the Kohn-Sham molecular electron density which, via the partitioning-specific weight function w i (r), is decomposed into its AIM densities {ρ i (r)}.

V i = r 3 ρ i (r)d 3 r ρ i (r) = w i (r)ρ(r) (2)
Once the set of static AIM polarizabilities in Eq.( 1) is obtained, a correspondent set of frequency-dependent ones is generated via Eq.( 3), where this time ω 0 j and C 0 6,j are the free atom characteristic excitation frequency and first dispersion coefficient.

α j (iν) = α j 1 -(iν/ω 0 j ) 2 ω 0 j = 4 3 C 0 6,j α 0 j 2 (3)
These frequency-dependent polarizabilities are, in the MBD@rsSCS model, gathered as diagonal elements of the frequency-dependent superpolarizability matrix A(iν) being one of the entries in the Dyson-like equation below which solution provides the screened superpolarizability matrix Ā(iν).

Ā(iν) = A(iν) -A(iν)T SR (iν) Ā(iν) (4)
The T SR represents a damped dipole-dipole interaction operator applied to the Coulombic interaction of two frequency-dependent spherical Gaussian charge distributions where its explicit expression, together with the one for A(iν), can be found in reference. [START_REF] Poier | N) stochastic evaluation of many-body van der Waals energies in large complex systems[END_REF] We note here that the Fermi damping function employed in the definition of T SR makes use of AIM van der Waals radii which can also be obtained by a volume rescaling similarly to what discussed for polarizabilities. [START_REF] Tkatchenko | Accurate molecular van der Waals interactions from ground-state electron density and free-atom reference data[END_REF] The solution of Eq.( 4) for a set of frequencies, and a consequent partial contraction of the converged { Ā(iν)}, gives a set of screened frequency-dependent atomic polarizabilities {ᾱ j (iν)} which are used to approximate the Casimir-Polder integral providing screened characteristic excitation frequencies {ω j }.

C6,j = 3 π ∞ 0 ᾱj (iν)ᾱ j (iν)dν ωj = 4 3 C6,j ᾱj (0) 2 (5)
The set of screened excitation frequencies as well as the screened static atomic polarizabilities define the MBD potential matrix shown in Eq.( 6) for a general ij block. T LR represents the range-separated damped dipole-dipole interaction matrix which explicit expression is also found in reference. 37

V ij = δ ij ω2 i + (1 -δ ij )ω i ωj ᾱi (0)ᾱ j (0)T LR ij (6)
The trace of √ V defines the interaction energy E int of the CFDs in the system 37 while its zero-point value E 0 is given by the sum of all screened excitation frequencies. Finally the difference between E int and E 0 gives the target MBD@rsSCS energy, Eq.( 7), which is coupled to the KS-DFT one to include non-additive dispersion contributions.

E MBD = E int -E 0 = 1 2 Tr[ √ V] - 3 2 N i=1 ωi (7)
In the original MBD@rsSCS model just briefly reviewed, E MBD is coupled to the molecular electron density via AIM volume partitioning introduced in Eq.( 2).

In this Letter instead we show that the explicit electron density partitioning can be avoided by learning AIM volumes via a DNN model without affecting the original MBD@rsSCS model's accuracy.

Bereau et al. and more recently, Mulhi et al. used ML on atomic volumes inside vdW model

to capture many body effects. [START_REF] Bereau | Non-covalent interactions across organic and biological subsets of chemical space: physics-based potentials parametrized from machine learning[END_REF][START_REF] Muhli | Machine learning force fields based on local parametrization of dispersion interactions: Application to the phase diagram of C 60[END_REF] Both have developed a Gaussian approximation potential (GAP) force field on TS polarizability rescaling. While GAP has shown to outperform neural networks in predicting energies with small-sized data set, e.g few thousands of data, its poor computational scaling O(N 3 ) prevents its use on very large training sets and thus to build a general purpose MBD model [START_REF] Zuo | Performance and cost assessment of machine learning interatomic potentials[END_REF] Finally, these models are either restricted to pairwise interactions or do not scale linearly with respect to the number of atoms as our Stochastic reformulation of the MBD equations was introduced only recently. [START_REF] Poier | N) stochastic evaluation of many-body van der Waals energies in large complex systems[END_REF] Isayev et al. [START_REF] Zubatyuk | Accurate and transferable multitask prediction of chemical properties with an atoms-in-molecules neural network[END_REF] recently extended their 5 million chemical conformations, the ANI-1 data set, with atomic volumes computed at the ωB97x/def2-TZVPP level with MBISA partitioning.

In virtue of its size and diversity, this data set is here employed in building our DNN to be coupled to the MBD@rsSCS model. Here we restrict ourselves to structure composed of only C, H, N and O, thus reducing the actual data set size to 4.6 million conformations.

In the MBISA weight function w i (r), each of the reference pro-atomic densities ρ 0 i (r) is expanded into m i Slater functions, m i being the number of shells of atom i placed at R i .

w i (r) = ρ 0 i (r) N j=1 ρ 0 j (r) ρ 0 i (r) = m i σ=1 N i,σ k 3 i,σ 8π exp - r -R i k i,σ (8)
In the scheme, the population N i,σ and width k i,σ of each shell are free-variables which are optimized so that the loss of information upon molecule formation is minimized. [START_REF] Verstraelen | Minimal basis iterative stockholder: atoms in molecules for force-field development[END_REF] To handle such large data set, a deep neural network is the natural choice. [START_REF] Westermayr | Long-range dispersion-inclusive machine learning potentials for structure search and optimization of hybrid organic-inorganic interfaces[END_REF] In particular, we use as machine learning model a feed-forward DNN with the ANI-like symmetry functions (SFs). [START_REF] Smith | ANI-1: an extensible neural network potential with DFT accuracy at force field computational cost[END_REF] The ANI's SFs are a subfamily of Behler-Parinello's ones [START_REF] Behler | Generalized neural-network representation of highdimensional potential-energy surfaces[END_REF] which traduce an atomic local environment i into an atomic environment vector (AEV)

G i = {G R i , G A i } where
G R i and G A i represent its radial and angular contributions respectively. Although SFs development is an intensive field of research and more accurate models have been developed since (ωACSF, 44 SOAP 45 among others), we stick to the ANI's original SFs as they were shown to successfully predict complex local properties such as, in the case of AIMNET, multipoles and volumes. [START_REF] Zubatyuk | Accurate and transferable multitask prediction of chemical properties with an atoms-in-molecules neural network[END_REF] Moreover, ANI's SFs have the great advantage of being computationally efficient as they rely on 2-body terms thus making the overall DNN model linear scaling with the system's size.

The DNN part of the combined DNN-MBD model relies on Scikit-learn, [START_REF] Pedregosa | Scikit-learn: machine learning in Python[END_REF] PyTorch 47 and TorchAni. [START_REF] Gao | TorchANI: A free and open source PyTorch-based deep learning implementation of the ANI neural network potentials[END_REF] They are all included in the Tinker-HP neural network module, which implementation will be detailed in a forthcoming dedicated paper (T. Jaffrelot Inizan et al., 2022).

We kept the original ANI's SFs parameters as we did not see major differences after tuning them. We empirically tested multiple neural network architectures (further details are found in the Supporting Information (SI) Figure 3) and the best performance was obtained with information. We used the Exponential Linear Units (ELU) activation function [START_REF] Clevert | Fast and accurate deep network learning by exponential linear units (ELUs)[END_REF] while the model's parameters were initialized with the so-called "He" initialization and updated with Hutter's AdamW algorithm during the training procedure. [START_REF] Loshchilov | Decoupled weight decay regularization[END_REF] Within the AdamW algorithm, the factor was set to 0.5 and the patience to 100. The initial learning rate was set to 10 -3

and the early stopping learning rate was set to 10 -6 . The ANI-1 dataset was shuffled and split into training and validation set containing 80% and 20% respectively of the full dataset.

The networks were trained for 6000 epochs with a batch size of 2560.

The ANI-1 data set, upon which our DNN model is trained, consists of AIM volumes computed at the ωB97x/def2-TZVPP level. The model is trained on volume ratios rather than pure AIM volumes as the narrower distribution of the former allows for a DNN's better performance without the need for re-scaling. Indeed, the atomic volumes ratio for C, H, O, N (see Figure 1 of the SI) is between 0.1 and 1. The DNN model providing AIM volumes' ratios is embedded in the Tinker-HP package where our linear-scaling and embarrassingly parallel stochastic MBD@rsSCS is also implemented. [START_REF] Poier | N) stochastic evaluation of many-body van der Waals energies in large complex systems[END_REF] The outcoming DNN-MBD model is coupled to the common semi-local PBE 52 functional as well as its hybrid PBE0 version 53 since this choices allow for comparisons with results ready available in literature. The optimal range-separation β parameters for both the PBE+DNN-MBD and PBE0+DNN-MBD methods are obtained by minimizing the mean absolute rel-ative error (MARE) on the widely employed S66x8 benchmark set consisting of 66 dimers placed at 8 different intermolecular distances for a total of 528 different structures where CCSD(T) interaction energies computed at CBS are used as reference.

All DFT computations employed Jensen's pcseg-3 basis set belonging to the family of segmented polarization-consistent 54 basis sets which, for DFT calculations, exhibits lower basis set errors than other gaussian basis sets as well as higher computational efficiency at given cardinal number as these basis sets were explicitly designed and optimized for DFT. [START_REF] Jensen | Unifying general and segmented contracted basis sets. Segmented polarization consistent basis sets[END_REF] Figure 2 shows the MARE as a function of the range separation parameter for PBE+DNN-MBD and PBE0+DNN-MBD methods.

The optimal β parameters are found to be 0.75 and 0.77 for the PBE+DNN-MBD and PBE0+DNN-MBD methods respectively. These values differ from the ones optimized for the original PBE/PBE0+MBD@rsSCS models [START_REF] Ambrosetti | Long-range correlation energy calculated from coupled atomic response functions[END_REF] and this has to be addressed to the different partitioning scheme employed. As pointed out by Vestraelen et al., [START_REF] Verstraelen | Minimal basis iterative stockholder: atoms in molecules for force-field development[END_REF] AIM densities computed via the Hirshfeld or HI partitioning, exhibit asymmetries i.e. they are aspherical with too much density in the bonding region. This density accumulation, relatively far away from the atomic nucleus, leads to larger values of radial moments, thus leading to larger AIM volumes compared to the ones obtained via the MBISA scheme (unaffected from this asymmetry artifact) for which less screening of volume-scaled AIM quantities (smaller β) is most likely to be needed. [START_REF] Verstraelen | Minimal basis iterative stockholder: atoms in molecules for force-field development[END_REF] We observe, nevertheless, that both PBE0+DNN-MBD and PBE0+MBD@rsSCS methods require a larger β parameter compared to their PBE corresponding models and this is consistent with the PBE0's improved description of short-range exchange-correlation effects due to the fraction of exact exchange included in the functional, as discussed in reference. [START_REF] Ambrosetti | Long-range correlation energy calculated from coupled atomic response functions[END_REF] The performance of the optimized PBE/PBE0+DNN-MBD methods is compared to different MBD models in terms of MAE and MARE for the S66x8 data set and the results are summarized in Figure 3 where actual values are reported in Table 1.

For the benchmark set here employed, the DNN-MBD model exhibits lower (although by a contained margin) errors both in its coupling to the PBE and PBE0 functionals compared to the standard MBD@rsSCS approach based on Hirshfeld AIM volumes as well the PBE+MBD@rsSCS/FI approach based on the fractionally ionic polarizabilities and HI AIM volume partitioning. For both the chosen functionals, the outcoming DNN-MBD model provides a mean absolute error in the S66x8 interaction energies which is below 0.25 kcal/mol compared to the reference CCSD(T) CBS golden standard.

To strengthen the analysis, we additionally computed the MAE and MARE for the S22 data set 57 by employing the range separation parameters previously optimized for the S66x8 set.

We can, in this way, employ the S22 set as a test set to validate our conclusions, Figure 3 (bottom) and Table 2.

Compared to the S66x8 set, the MAE and MARE values of our proposed PBE/PBE0+DNN-MBD models are, for the S22 set, higher however this is not surprising as no β optimization was performed this time. Let's note that all methods present errors that are larger in the Having been trained on a large and diverse set of AIM volumes, the outcoming DNN-MBD model inherits the strengths of the MBISA scheme discussed earlier in this Letter while completely bypassing the explicit density partitioning with a consequent decrease of the computational cost. We also note that the DNN model could be successfully trained with different AIM partitioning schemes due to the locality of the target quantities (volumes).

The presented density-free DNN-SMBD model is included in the Tinker-HP package [START_REF] Lagardère | a massively parallel molecular dynamics package for multiscale simulations of large complex systems with advanced point dipole polarizable force fields[END_REF] and will be released with the next version of the software. There, it can benefit from the linearscaling embarrassingly parallel performances of our stochastic formulation (SMBD) of the MBD key equations which remarkable computational performances have been recently discussed. [START_REF] Poier | N) stochastic evaluation of many-body van der Waals energies in large complex systems[END_REF] We believe that the present DNN-SMBD model can be beneficial in applications of dispersioncorrected DFT to large complex systems requiring an accurate yet extremely efficient inclusion of MBD effects. The DNN model, by avoiding the direct solution of the KS equations due to its electron density-free features, allows for the ready application of DNN-SMBD approach in the development of accurate ab initio-based force fields [START_REF] Gresh | Anisotropic, polarizable molecular mechanics studies of inter-, intra-molecular interactions, and ligandmacromolecule complexes. A bottom-up strategy[END_REF][START_REF] Naseem-Khan | Development of the Quantum Inspired SIBFA Many-Body Polarizable Force Field: Enabling Condensed Phase Molecular Dynamics Simulations[END_REF] and neural networks methodologies.

6 .

 6 Free atom volumes are computed at the same level as AIM ones. The correlation plots between the DNN model and the ab initio validation set reference is depicted in Figure 1. The root-mean-square-error (RMSE) and mean-absolute-error (MAE) are respectively 0.012 and 0.008 which is much less than the smallest value of the data set showing the good accuracy of our model. The final DNN model and the dataset used for the training can be download directly via the Zenodo repository located at the following address. 51

Figure 1 :

 1 Figure 1: Atomic volume correlation plot comparing the DNN prediction to DFT reference calculations for 1/100 of the validation set. The color bar scale reflects the density of points and correlate with the atomic volumes ratio distribution (Figure 2 of the SI).

Figure 2 :

 2 Figure 2: MARE (%) as a function of the range separation parameter for the PBE+DNN-MBD and PBE0+DNN-MBD methods.

Figure 3 :

 3 Figure 3: MARE (%) and MAE (kcal/mol) of PBE, PBE0, PBE+D3 56 and different MBD models (MBD@rsSCS, 21 MBD@rsSCS/FI 30 ) including our DNN-MBD for the S66x8 (top) and S22 (bottom) data sets.

Table 1 :

 1 [START_REF] Goerigk | Benchmarking Density Functional Methods against the S66 and S66x8 Datasets for Non-Covalent Interactions[END_REF] and different MBD models (MBD@rsSCS,[START_REF] Ambrosetti | Long-range correlation energy calculated from coupled atomic response functions[END_REF] MBD@rsSCS/FI 30 ) including our DNN-MBD for the S66x8 (top) and S22 (bottom) data sets. case of the S22 set compared to S66x8 (see Table1 and Table 2). Indeed, there are reasons for that and we can stress that the dimers employed in the S22 set are placed at equilibrium while the S66x8 set includes out of equilibrium dimers. In our case, the DNN-MBD model trained on S66x8, appears less biased towards equilibrium structures. Overall, as one can see from Table2, our DNN-MBD remains highly transferable and, with an error below 0.43 MAE (kcal/mol) and MARE(%) relative to the S66x8 data set for our DNN-based models as well as for few other dispersion correction ones. For the MBD-based models, the method-specific range separation parameter reported in parentheses refers to the one optimised for the S66x6 set. MAE and MARE are computed taking revised CCSD(T) CBS energies.

	kcal/mol compared to the reference CCSD(T) CBS golden standard, outperforms previous
	S22 results obtained with others methods.
	Model	MAE[kcal/mol] MARE%
	PBE	1.55	65
	PBE0	1.48	65
	PBE+D3	0.44	n.a.
	PBE+MBD@rsSCS (β = 0.83)	0.32	10.6
	PBE0+MBD@rsSCS (β = 0.85)	0.30	9.2
	PBE+MBD@rsSCS/FI (β = 0.83)	0.28	9.0
	PBE+DNN-MBD (β = 0.75)	0.25	9.0
	PBE0+DNN-MBD (β = 0.77)	0.23	6.9

Table 2 :

 2 MAE (kcal/mol) and MARE(%) relative to the S22 data set for our DNN-based models as well as for few other dispersion correction ones. For the MBD-based models, the method-specific range separation parameter reported in parentheses refers to the one optimised for the S66x6 set. MAE and MARE are computed taking revised S22 energies where, compared to the original paper, a larger basis set was employed.[START_REF] Takatani | Basis set consistent revision of the S22 test set of noncovalent interaction energies[END_REF]

	Model	MAE[kcal/mol] MARE%
	PBE	2.66	58
	PBE0	2.44	55
	PBE0+MBD@rsSCS (β = 0.85)	0.55	8.5
	PBE+MBD@rsSCS (β = 0.83)	0.49	8.9
	PBE+D3	0.48	n.a.
	PBE0+DNN-MBD (β = 0.77)	0.43	5.6
	PBE+DNN-MBD (β = 0.75)	0.41	6.6

Acknowledgement

This work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 810367), project EMC2 (JPP). Computations have been performed at GENCI (IDRIS, Orsay, France and TGCC, Bruyères le Chatel) on grant no A0070707671.

Supporting Information Available


M R Pereira
email: cristina.todisco@tetratech.com

FABRIC CHARACTERISATION IN TRANSITIONAL SOILS

Keywords: fabric, MIP, statistical parameters, transitional soils

A "transitional" mode of soil behaviour implies that dense and loose samples do not converge towards the same volumes within the strains and stresses applied by simple oedometer and triaxial tests. As this behaviour involves soils with different gradings and mineralogies (e.g. gap graded, well graded and/or mixed mineralogies), identifying the factors responsible is difficult. Nevertheless, it has been previously speculated that strong forms of fabric that are difficult to break down as strains and stresses are applied, might be the common cause.

This paper aims at investigating some elements of fabric at the microscale of transitional soils.

A gap graded and two well graded mixtures with large amounts of non-plastic fines were investigated by oedometer and triaxial tests. As it would be difficult to identify experimentally many commonly used elements of fabric in these soils, e.g. the contact network, mercury intrusion porosimetry (MIP) was used as a first step to characterise the evolution of pore size distributions (PSDs) of dense and loose samples undergoing the same stress paths, using the PSDs as a proxy of fabric. Multi-directional bender element testing was performed to confirm the isotropy of the elastic stiffness, from which it might be inferred that the fabric is also isotropic. Statistical parameters of the PSDs were calculated, the changes of which were related to the evolution of macroscale void ratios.

The robust fabrics causing lack of convergence were characterised by a complex evolution of the PSDs, the initial differences of which could not be erased during conventional testing. This work also provided a simple method to examine the fabric of particularly well graded or gap graded materials, for which other techniques, such as CT or SEM, could not reveal the multiscale nature of the fabric.

INTRODUCTION

A range of soils has now been observed to have a so called "transitional" mode of behaviour, for which convergence of loose and dense samples towards unique volumes is not seen either in compression or shearing within the range of strains that may be applied by simple oedometer or triaxial tests. The factors responsible for this have not clearly been identified but it has been speculated that it results from strong fabrics at the microscale that are difficult to break down [START_REF] Coop | Limitations of a Critical State framework applied to the behaviour of natural and "transitional" soils[END_REF]. However, identifying those elements of fabric responsible has proven elusive, mainly because of the difficulty of defining the fabrics of soils composed of a wide range of particle sizes and/or different mineralogies and that may have undergone complex geological processes (e.g. [START_REF] Ferreira | Problems in identifying the effects of structure and critical state in a soil with a transitional behaviour[END_REF][START_REF] Altuhafi | The mechanics of subglacial sediment: an example of new "transitional behaviour[END_REF]). Some natural soils characterised by strong forms of fabrics may show analogous behaviour of transitional soils. These robust fabrics can be observed at different scales and in different forms and due to this variability they are often classified as being heterogeneous. Heterogeneity might be found in particle and pore arrangements and topology and in force chain transmission, although the latter has been less investigated experimentally. For example, natural alluvial clayey soils (e.g. [START_REF] Coop | The compression of sediments at the archaeological site of Sibari[END_REF][START_REF] Cotecchia | One-dimensional compression of a natural clay: structural changes and mechanical effects[END_REF]) often have heterogeneous fabrics at the mesoscale with thin depositional layers of fine and coarse soils and have been found to have compression behaviour that is not convergent with the compression lines of their remoulded soils. Other clays may have silt-sized aggregates formed by clays minerals and are heterogeneous at the microscale (e.g. [START_REF] Fearon | Reconstitution: what makes an appropriate reference material?[END_REF][START_REF] Atkinson | Destructuring and disaggregation of Mercia Mudstone during full-face tunnelling[END_REF]), with different behaviours according to the degree of aggregate destructuration applied. Also DEM simulations on fractally graded sand mixtures showed transitional mode of behaviour when subjected to conventional laboratory compressive stress levels (<8MPa) [START_REF] Minh | A DEM investigation of the effect of particle-size distribution on one-dimensional compression[END_REF].

In this case robust fabrics were observed by analysing the force chain transmission. Strong force chains were carried by big particles, with the large voids compressed, while the small particles were either weakly loaded by the adjacent big particles or filling voids without transmitting any force, so that the small voids were little affected.

Based on these observations and the fact that their soils had isotropic strains, Shipton and Coop [START_REF] Shipton | Transitional behaviour in sands with plastic and non-plastic fines[END_REF] speculated that the fabric that was responsible for transitional behaviour might have a heterogeneous rather than anisotropic nature at the microscale, although the representative element volume for that micro fabric could not be defined. It is possible that the former is a more robust characteristic, since the latter may often be erased or at least modified as the test proceeds. The effects of fabric anisotropy on soil behaviour have been extensively studied at the macro and mesoscales [START_REF] Oda | Initial fabrics and their relations to mechanical properties of granular material[END_REF][START_REF] Kuwano | On the applicability of cross-anisotropic elasticity to granular materials at very small strains[END_REF][START_REF] Yang | Undrained anisotropy and rotational shear in granular soil[END_REF], but there is very much less research on the effects of fabric heterogeneity at the microscale. X-ray CT scanning and scanning electron microscopy, SEM, have both been used extensively to characterise soil structure and relate changes to the macromechanical behaviour (e.g. [START_REF] Delage | Study on the structure of a sensitive Champlain clay and of its evolution during consolidation[END_REF][START_REF] Carraro | Shear strength and stiffness of sands containing plastic and nonplastic fines[END_REF][START_REF] Hattab | Experimental study of kaolin particle orientation mechanism[END_REF][START_REF] Fonseca | Quantifying the evolution of soil fabric during shearing using directional parameters[END_REF][START_REF] Fonseca | Quantifying the evolution of soil fabric during shearing using scalar parameters[END_REF]). But for transitional soils, which are often gap graded or very well graded, a clear detection of the fabric elements responsible is difficult to achieve due to the difficulty of examining the fabrics at the scales of the smaller and larger particles simultaneously (see Appendix Figure A1 for an example). Nocilla et al. [START_REF] Nocilla | The mechanics of an Italian silt; an example of "transitional" behaviour[END_REF] and Shipton and Coop [START_REF] Shipton | Transitional behaviour in sands with plastic and non-plastic fines[END_REF] both tried SEM in unsuccessful attempts to identify the fabric responsible for the transitional behaviour that they observed. Mercury intrusion porosimetry MIP overcomes this shortcoming, detecting pore sizes from a few nm to a few hundreds of μm, although it is limited to characterising the soil structure only in terms of pore size distribution (PSD).

This paper aims at investigating the fabrics of transitional soils by means of MIP, examining the PSD and its evolution during conventional laboratory testing. In MIP tests, the volume of voids is measured over the volume of sub-samples of about 0.5-1cm 3 . It was found that the Representative Elementary Volume REV for an unsaturated fine sand was about 30-45 times larger than the particle size, equal to 10-15mm [START_REF] Bruchon | X-ray microtomography characterisation of the changes in statistical homogeneity of an unsaturated sand during imbibition[END_REF]. The REV for the materials presented in this work was not investigated but given their grain size distributions and the reasonable repeatability of most of the MIP results, it was assumed that the size of the MIP samples was large enough to capture representative PSDs of the whole samples. The PSDs were analysed by statistical parameters that were related to the macromechanical behaviour of the soils. It should be emphasised that MIP can only give information about one aspect of soil fabric, i.e. the pore size distribution, and gives no information about the precise nature of particle and void distributions and orientations that create the fabric. Since MIP cannot address fabric anisotropy, this has been inferred indirectly by means of multi-directional bender elements.

MATERIALS AND METHODS

The effects of the fabric have been studied for three mixtures, of which the mechanical behaviours were described in detail in Todisco and Coop [START_REF] Todisco | Quantifying rates of volumetric convergence in soils with complex gradings[END_REF] and Shipton and Coop [START_REF] Shipton | On the compression behaviour of reconstituted soils[END_REF][START_REF] Shipton | Transitional behaviour in sands with plastic and non-plastic fines[END_REF], and additional tests were carried out specifically to examine the fabric. The SPF (sand with plastic fines) is a gap graded mixture made of 75% Thames Valley sand [START_REF] Takahashi | Assessment of standard research sand for laboratory testing[END_REF] and 25% kaolin.

Its mechanical behaviour revealed apparently parallel normal compression lines (NCLs) and critical state lines (CSLs) in the state plane for different initial void ratios [START_REF] Shipton | On the compression behaviour of reconstituted soils[END_REF][START_REF] Shipton | Transitional behaviour in sands with plastic and non-plastic fines[END_REF]. The LBS and LMS were very well graded samples of Leighton Buzzard quartz sand and a crushed limestone from China. In both cases the maximum particle size was 600 μm, and within the sand fraction a fractal grading of 2.57 was used [START_REF] Altuhafi | Changes to particle characteristics associated with the compression of sands[END_REF][START_REF] Todisco | Quantifying rates of volumetric convergence in soils with complex gradings[END_REF]. Since it was not possible to control the grading within the fines fraction also to be fractal, the gradings were completed with 40% of crushed quartz silt or crushed limestone silt. The grain size distributions of the mixtures are shown in Fig. 1. The SPF samples were created by the moist-tamping method, varying the initial water content and number of layers as indicated in Shipton and Coop [START_REF] Shipton | Transitional behaviour in sands with plastic and non-plastic fines[END_REF] to change the initial void ratio. The samples of LBS and LMS were created by the dry compaction method which follows the procedure of the under-compaction method [START_REF] Ladd | Preparing test specimens using undercompaction[END_REF] but using dry reconstituted soils.

The one-dimensional tests on SPF were performed in conventional oedometers, using a 50mm diameter ring to reach a vertical stress level of about 8MPa. Smaller 20 and 30mm rings were used to reach vertical stresses of around 20-50MPa, but these had a floating ring design to minimise the side friction. The triaxial tests were carried out in typical stress path type apparatuses. The samples were saturated by first circulating CO2, then flushing with de-aired water and finally by applying back pressures of at least 200kPa, obtaining B values of 0.96-0.98. After connecting a suction cup to link the axial loading system rigidly to the sample [START_REF] Atkinson | Discussion on: The measurement of soil stiffness in the triaxial apparatus[END_REF], isotropic compression to different stress levels was followed by drained shearing under axial strain control, typically using a gradual increasing rate from 0.05%/h in the small strain region (<0.1% axial strain) to 0.4%/h at large strains. This was a pragmatic choice to complete the tests in a relatively short time, while ensuring complete drainage, although very small rate effects in the stress-strain curves were observed. However, the samples were retrieved at the end of the tests after reaching the critical state. It has been shown that rate effects become negligible on both stress and state planes as axial strains increase (e.g. [START_REF] Tatsuoka | Time-dependent shear deformation characteristics of geomaterials and their simulation[END_REF][START_REF] Li | Strain rate dependence of the critical state line of reconstituted clays[END_REF]). Full details of the tests are given in Tables 1 and2.

Multi-directional bender element testing

A Bishop and Wesley [START_REF] Bishop | A hydraulic triaxial apparatus for controlled stress path testing[END_REF] triaxial apparatus was fitted with T-configuration lateral bender elements [START_REF] Pennington | Anisotropy of G0 shear stiffness in Gault clay[END_REF] able to measure the stiffnesses Ghv and Ghh. These were inserted through the membrane using a specially designed mould, described in detail in Todisco [START_REF] Todisco | Behaviour of non-convergent soils: effect of particle size, mineralogy and fabric[END_REF]. The pedestal and top platen also housed axially orientated bender elements to measure Gvh. The data from these were consistent with Ghv and Ghh, but since the vertical bender elements have different boundary conditions to the lateral, data from them tended to increase the scatter and so for this reason they have not been presented in the analysis. The shear wave velocities were calculated using the first arrival time [START_REF] Abbiss | Shear wave measurements of the elasticity of the ground[END_REF] ensuring a consistent choice over a range of frequencies from 8 to 20 kHz.

Mercury intrusion porosimetry (MIP) tests

The MIP tests were carried out using an AutoPore IV 9500. Applying Eq. 1 [32] it was possible to obtain the pore diameter intruded by the pressurised mercury 𝑝 = -

1
where p is the absolute pressure applied to the mercury, n is a coefficient accounting for the pore shape, a value equal to 4 (corresponding to a cylinder) being adopted in this study, σHg is the surface tension of the mercury equal to 0.484 N/m at 25° C, θ is the contact angle between the pore contour and the mercury and in this study is assumed equal to 130° and d is the pore throat diameter. Different contact angles were reported in the literature with values ranging from 130° [33] to 160° [START_REF] Diamond | Pore size distribution in clays[END_REF][START_REF] Penumadu | Compressibility effect in evaluating the pore-size distribution of kaolin clay using mercury intrusion porosimetry[END_REF]. The SPF samples were freeze-dried after being carefully waxcoated at the end of the tests, while the LBS and LMS samples were oven-dried. The freezedrying technique was used for the sand-clay mixture SPF to avoid bulk shrinkage and changes in pore size distribution [START_REF] Diamond | Pore size distribution in clays[END_REF][START_REF] Delage | Study on the structure of a sensitive Champlain clay and of its evolution during consolidation[END_REF], while it was not necessary for the LBS and LMS because they were non-plastic.

In each case the MIP samples were carefully trimmed to a size roughly equal to 1cm 3 , but optimised depending on the porosity of each sample to obtain the best resolution from the apparatus. The oedometer and triaxial sub-samples were retrieved after unloading. The trimming of the SPF was relatively straightforward, but the trimming and handling of the sand samples required extreme care and neither sand could be sampled and tested in its initial state, the samples either collapsing during trimming or during initial immersion in the mercury. For the final states, some, but not all tests were successful. The successful tests resulted perhaps because their very well graded nature ensured that they just had sufficient particle interlock after loading to allow the test to be done successfully, or perhaps because there was a tiny amount of bonding created by loading, especially in the LMS. The slight cohesion in the samples could not have resulted from suction as the MIP test is carried out under a very high vacuum.

The pore size distributions, PSDs, are shown with the x-axis on a logarithmic scale, as often adopted for data spreading over many orders of magnitude. PSDs were analysed in terms of statistical parameters that offered a valuable tool to describe the shape of a probability density function. Sedimentologists use this approach to describe particle size distributions by reading selected percentiles of the cumulative curves (e.g. [START_REF] Folk | Brazos River bar: a study in the significance of grain size parameters[END_REF]). Here, the method of the scaled moments [START_REF] Pearson | Editorial note to "Inequalities for moments of frequency functions and for various statistical constants[END_REF] was preferred as it can be applied to the majority of the curve shapes (i.e. nonnormal distributions). The statistical parameters of the PSDs were the mean (μ), standard deviation (σ), skewness (γ) and kurtosis (κ). An i-th moment is scaled when it is divided by the standard deviation to the i-th exponent. The first scaled moment is 0 because the moment with exponent 1 around the mean is 0, the second is 1 because the moment with exponent 2 around the mean is the variance σ 2 , the third and the fourth scaled moments are skewness and kurtosis.

The statistical parameters were calculated manually using the discrete values of the functions obtaining from the MIP tests. Equations 2-6 explain the procedure in detail. At each discrete value of xi that is the log of pore diameter, the i-th area under the probability function ΔAi is equal to:

∆𝐴 = 𝑓(𝑥 )∆𝑥 2
where f(xi) is the i-th value of the density function, also called incremental pore volume in the following figures. The mean was calculated as

𝜇 = ∑ ∆ ∆ ∑ ∆ 3
where Δxi-0 is the distance of the i-th log pore diameter interval from the origin. The standard deviation σ was calculated as:

𝜎 = ∑ ∆ (∆ ) ∑ ∆ 4
where Δxi-μ is the distance of the i-th log pore diameter interval from the mean. The values shown in the Results section are the inverse of the logarithmic values obtained by Eq.4. The skewness γ and the kurtosis κ were calculated as:

𝛾 = ∑ ∆ (∆ ) ∑ ∆ 5 𝜅 = ∑ ∆ (∆ ) ∑ ∆ 6
The statistical parameters may be compared to those of a normal distribution that has a γ of 0 and κ of 3. The value of γ locates the centre of mass of the distribution, a negative value defining a left-skewed distribution with the centre of mass to the left of the mean and longer tail towards the right. The value of κ defines the sharpness of the peak and the thickness of the tails; if it is greater than 3 then the peak is sharper and tails longer and thicker than for a normal distribution.

These statistical parameters for the PSDs were plotted against the final void ratios of the tests to try to relate the macromechanical behaviour to the fabric.

RESULTS

Changes to void ratio in compression and shear

Figure 2 shows the mechanical behaviour of selected SPF, LBS and LMS samples that were subjected to MIP and bender elements (BE) testing. These were part of a more extensive experimental campaign, which investigated the transitional behaviour of these mixtures in compression and shearing [START_REF] Todisco | Quantifying rates of volumetric convergence in soils with complex gradings[END_REF] although additional tests have been carried out in this investigation of fabric. The sample names help to indicate the initial void ratio and the maximum stress level. For example, LI and DI indicate loose and dense samples of SPF in their "initial state" which was one of one-dimensional compression to about 50kPa, so that the samples were firm enough to be handled and trimmed; LF and DF were different samples with initial void ratios similar to LI and DI but compressed to 8MPa and then retrieved after unloading to 50kPa for the MIP tests.

Values of mean effective stress p' were plotted for the oedometer tests on SPF and LMS by assuming k0=1-sinφ' [START_REF] Jaky | Pressure in silos[END_REF], where k0 is the coefficient of earth pressure at rest equal to σ'3/ σ'1

for zero lateral strain and φ' is the angle of shearing resistance. In standard oedometer tests only the vertical stress σ'1 is known, calculated from the applied load, while the horizontal stress σ'3 is obtained by multiplying σ'1 by k0. The mean effective stress p' is the first stress invariant equal to , where σ'1 and σ'3 (= σ'2 i.e. axisymmetric conditions) are the major and minor principal effective stresses. For the LBS two samples were tested using lubricated end platens in the triaxial, but as discussed by Todisco and Coop [START_REF] Todisco | Quantifying rates of volumetric convergence in soils with complex gradings[END_REF] this did not affect the data significantly in comparison to the much larger differences of void ratio between different samples. The samples tested using bender elements (BE) and MIP are indicated in the graphs.

Multi-directional bender element tests

First, the characterisation focused on the fabric anisotropy and the data from the multidirectional bender element tests during isotropic compression are given in Fig. 3. These indicated that the elastic stiffnesses were isotropic, the Ghv and Ghh values being very similar, from which it is inferred that the fabrics were also isotropic. The BE tests on SPF confirm the isotropy that was suggested by Shipton and Coop [START_REF] Shipton | Transitional behaviour in sands with plastic and non-plastic fines[END_REF] from an examination of the axial and volumetric strain increments during isotropic compression. A loose (LBS16) and dense (LBS17) sample of LBS shows little difference in the stiffnesses. From small strain probes Shipton and Coop [START_REF] Shipton | Transitional behaviour in sands with plastic and non-plastic fines[END_REF] had tentatively reached a similar conclusion for the SPF, but their data were rather more scattered.

MIP results

The density distribution curves for the tests on SPF are given in Figure 4. The samples in their initial states (post-50kPa-compression), LI and DI show that they have significantly different PSDs, the loose (LI) having larger pores and a broader peak than the dense. The broad peak in fact consists of a slight double peak. The small trough at about 6μm may be disregarded as it occurs at the transition from the low to the high pressure analysis port. Compression moved the PSDs towards the region of smaller pores, as expected but it did not influence the distributions of pores smaller than 0.2m. The broad peak of loose sample LI was reduced in the final state (LF) as the individual peak for the larger pores reduced while that for the smaller pores increased. For the final dense sample (DF) there is a very much increased difference between the two initial peaks (DI), with the peak for the larger pores reduced very much more than for the loose samples. The initial dense distribution (DI) and final loose (LF) are actually quite similar, and on Fig. 2 they have quite similar void ratios. The final PSDs at the same vertical stress of about 8MPa remain significantly different, corresponding to their different void ratios, the denser sample having far fewer large voids. As the PSD is a fabric element, it could be concluded that these fabrics, as characterised here solely by the pore distributions, are robust and cannot be erased completely by compression. None of the PSDs showed the very marked bimodality found by Juang and Holtz [START_REF] Juang | Fabric, pore size distribution, and permeability of sandy soils[END_REF], who tested a similar mixture of 30% kaolin-70% Ottawa sand. Perhaps, the compaction method of Juang and Holtz [START_REF] Juang | Fabric, pore size distribution, and permeability of sandy soils[END_REF] generated a different fabric to the moist tamping method used here.

Figure 4b shows the PSDs of final samples of SPF sheared drained in the triaxial at 300kPa. It is possible that the data might be affected by experimental uncertainties, exacerbated by the small number of tests. However preliminary conclusions can be drawn, which seem consistent with the mechanical behaviour, although a more exhaustive validation is needed from future research. The distributions are quite different to those for the oedometer tests, which might be because of the different strains applied during one-dimensional compression and shearing and to the different volumetric behaviours of the loose and dense samples; SPF1 was contractive while SPF2 was dilative. Nevertheless, the distributions for the loose (SPF1) and dense samples (SPF2) are again quite different, which may justify their differences in void ratio in Fig. 2a.

The fabrics, in terms of PSDs, have again not converged even after shearing to about 30% axial strain.

To check whether fabric can be related to convergence, two control tests were carried out, testing dense and loose kaolin samples in the oedometer. These reached a unique NCL at about 100kPa (Fig. 5a). The samples were prepared as slurries at different initial water contents in order to obtain different initial void ratios. In Fig. 5b, the PSDs soon after compression have peaks at almost the same pore diameter (about 0.2μm) and are very similar over the full range of pore sizes.

Because of the difficulties in trimming and carrying out the tests on the sands, many of the MIP tests were not successful. A selection of those that were comparable is shown in Fig. 6 where the stress levels and final void ratios have been indicated for each sample. All the PSDs can be defined as unimodal independently of the loading type. Unfortunately, the unimodal shapes highlight that the differences for the sands are much less clear than for SPF only showing that dense samples have smaller modes (most frequent value of pore diameter) and extra smaller pores than the loose ones, the PSDs being slightly shifted to the left. For the LMS, the two oedometers were loaded to the same vertical stress of about 50MPa and the compression paths of the loose (LMS-OED1) and dense (LMS-OED2) samples were tending to converge slowly (Fig. 2c) but there were still significant differences in both void ratio and PSDs. The PSDs of the oedometers resemble a normal distribution, except for some lack of symmetry of the tails due to the presence of large pores.

The PSDs of the mixtures, especially the SPF ones, show that compression to stress levels smaller than 8MPa affects only pore between 0.2 and 10m, leaving unchanged the distributions of the smaller ones. Although MIP tests cannot investigate force chain transmission, it might be inferred that transitional behaviour in the mixtures arises because forces are not carried homogeneously by all the particles [START_REF] Minh | A DEM investigation of the effect of particle-size distribution on one-dimensional compression[END_REF]. This justifies the changes of PSDs only in specific regions of pores.

Overall, the results of MIP tests tend to confirm that different PSDs are associated with different void ratios for the type of transitional soils presented in this work, while MIP tests on convergent samples of kaolin reached a unique PSD.

Statistical analysis of the PSDs

In Fig. 7, the mean values (solid markers) of the PSDs of the oedometers and triaxials on SPF are rather different but both decrease as void ratio decreases. It seems that for the oedometers, where there are more data, the mean is fairly well related to void ratio, no matter whether it is at the start or end of the test. The direct comparison between the final values of LF and DF,

shows that the final mean does not converge to a unique value. In contrast the standard deviation of the oedometers does not vary significantly, and is similar for the initial and final values and for dense and loose samples. However, that of the triaxials decreases with decreasing void ratio. The standard deviation is larger than the mean probably because the pore size covers several orders of magnitude. The skewness γ and kurtosis κ of the oedometer and triaxial samples are more similar, but they are both slightly lower for the triaxials. With some data scatter, it again seems to be the case that for the oedometers the skewness and kurtosis are related to void ratio but both increase with decreasing void ratio, so again the final dense and loose samples have distinctly different values. In summary, after compression and shearing the PSDs of the SPF evolve into those with smaller pores on average, longer and thicker tails in the region of small pores >0.2m with a sharper peak (increasing kurtosis, κ) and a centre of mass located in the region of large pores (right skewed, increasing positive γ).

Figure 8 shows the statistical parameters for LBS and LMS. The differences are not large, mostly because the sample void ratio differences were also much smaller than for the SPF and the data are few. The mean of both sands decreases as void ratio decreases, as might be expected, but much less than for the SPF. This might be because the sand-clay mixture is more compressible overall than the well graded sands. In this case, the standard deviations increase slightly and the skewness of LBS samples increases as void ratio reduces, like the SPF but here κ remains nearly constant. The constant value of κ of LBS indicates that the sharpness of the peaks and the thickness of the tails are not significantly different between the dense and loose final samples.

Both the γ and κ of LMS decrease significantly in contrast to what was seen in the other soils.

The mean and skewness do not change as much as observed in the SPF samples, even allowing for the smaller differences on void ratio.

Table 3 summarises the statistical parameters of comparable samples of the various soils, each comparison with a similar type of test and stress level. The median and the diameter ratio d60/d10 considering the pore sizes at the 60 and 10% percentiles of the cumulative distributions have been added as a quantification of the symmetry and uniformity of the data. The ratio d60/d10 is similar the coefficient of uniformity adopted for the characterization of grain size distributions.

If it is less than 4, the soil is defined as uniformly graded. Generally, the mean values for the gap graded SPF vary much more between loose and dense samples than those of the sands, but this is expected since the SPF void ratios cover a wider range, both at the start and end of tests.

The median decreases with decreasing void ratio and shows generally smaller values than the mean. This indicates that the distributions are not symmetric but shifted toward the region of large pores, i.e. positive skewness values. The standard deviation of the SPF is larger than that of the sands, so that the SPF has more variability around the mean diameter than the well graded sands. The skewness is always positive and does not vary greatly between the well graded sands and the sand with plastic fines, the values ranging between 0.63 and 1.88. The kurtosis values are generally larger than 3, with the exception of test SPF1. This means that all the PSDs have sharper peaks and thicker tails than a normal distribution. The ratio d60/d10 is larger than 4 for SPF and LBS samples indicating that the pores are poorly sorted while LMS samples show more uniform distributions. As the samples become denser the uniformity of the distributions increases revealing that the compression and/or shearing tend to reduce the large initial differences in pore size.

CONCLUSIONS

The fabrics of three soils that might be described as "transitional" were characterised in terms of their pore size distributions, which were characterised as heterogeneous because covering large range of pore sizes. The similarities between stiffnesses, Ghv and Ghh, measured by multidirectional bender elements, suggested that in each case the fabrics are isotropic. The stiffnesses of the LBS were found to be poorly related to the densities of the samples and mostly dependent on the stresses applied. These results emphasise that the transitional mode of behaviour found for the three mixtures, or lack of convergence of void ratios in simple laboratory tests, must be a real soil behaviour, and is not related to any inherent anisotropy that might be generated during sample preparation.

MIP testing proved to be a good technique to characterise transitional soil behaviour in terms of micro fabric changes (PSDs changes). The MIP tests on SPF showed that the robust initial PSDs of the dense and loose samples were distinctly different and not erased during compression and/or shearing. Also pores between 0.2 and 10μm experienced the largest changes during compression, while the smaller ones remained almost unaffected. The MIP tests on the LBS and LMS were very much more difficult to conduct, with many failed tests and more scattered data. Also the narrower range of void ratios achieved, the unimodal nature of the PSDs and smaller changes of PSD with void ratio in the case of LBS, meant that they were more difficult to characterise than the SPF. Nevertheless, the data allow some significant conclusions still to be drawn. The loose and dense LBS samples sheared at the same isotropic pressure and also the oedometer tests on LMS showed different final PSDs that corresponded to their different final void ratios, the denser samples having smaller pores than the loose ones, as expected, but also with other changes to the shapes of the distributions. Two control oedometer tests on samples of kaolin, confirmed that a convergent void ratio can correspond to a convergent PSD. This might indicate that void ratios and pore distributions are directly linked for these soils: different void ratios correspond to different PSDs and similar void ratios have similar PSDs. However, in other cases it is not excluded that soils might have similar void ratios but different PSDs (e.g. [START_REF] Juang | Fabric, pore size distribution, and permeability of sandy soils[END_REF]).

Statistical parameters were applied to the PSDs for the first time. The quantitative parameters were linked to the state of the soils, providing a description of the evolution of the PSDs during compression and/or shearing. The mean and median of all the mixtures decreased as the samples became denser, while the trends of standard deviation, skewness, kurtosis and d60/d10 depended on the type of soil and test. This work offers a new insight on the fabric of transitional soils. The robust fabrics causing lack of convergence were isotropic and characterised by a complex evolution of the pore size distributions, the initial differences of which could not be erased during conventional testing. It also provided a simple method to examine the fabric of particularly well graded or gap graded materials, for which other techniques, such as CT or SEM could not reveal the multi-scale nature of the fabric. The statistical analysis of the PSDs and their relationship to soil states might be adopted more widely for other materials.

Furthermore, such analyses open new perspectives on modelling the behaviour of soils based on the knowledge of their PSD and its evolution with state parameters. For instance, the present work lets envisage an extension to soils of the work by Arson and Pereira [START_REF] Arson | Influence of damage on pore size distribution and permeability of rocks[END_REF] and Pereira and

Arson [START_REF] Pereira | Retention and permeability properties of damaged porous rocks[END_REF], who related the hydro-mechanical behaviour of damaged rocks to their PSD.

TABLES

Figure 1 Figure 2 Figure 3

 123 Figure 1 Grain size distributions of the tested mixtures

Figure 4 Figure 5 Figure 6 Figure 7 5 skewnessFigure 8

 456758 Figure 4 Density functions of the intruded volume of mercury of a) oedometer and b) triaxial samples of SPF.

Figure

  Figure A1 CT scan images of oedometer tests on LBS compressed to 50kPa: a) loose and b) dense (beam energy 21keV, average voxel resolution 0.625 μm).

Table 1

 1 Details of oedometer tests on SPF and LMS

	Test	Initial	Void ratio at	Final void	Vertical	p'
		void ratio	maximum σ'v	ratio	[kPa] stress, σ'v	[kPa]**
	LI	0.893	0.501	0.501	56	37 (k0=0.49)
	DI	0.429	0.371	0.371	56	37 (k0=0.49)
	LF	1.068	0.306	0.357	7700	5000 (k0=0.49)
	DF	0.443	0.263	0.306	7700	5000 (k0=0.49)
	Loose-kaolin	2.050	0.512	0.788	7700	-
	Dense-kaolin	1.804	0.512	0.747	7700	-
	LMS-OED1*	0.569	0.246	0.256	48000	28000 (k0=0.38)
	LMS-OED2*	0.420	0.152	0.196	48000	28000 (k0=0.38)

*selected samples for MIP tests, ** k0 values calculated from Jaky's relation (k0=1-sinφ')

Table 2

 2 Details of triaxial tests on SPF, LBS and LMS

	Test	Initial	Void ratio	Void	p' end of	p' end of	Volumetric	Shear
		void ratio	end of iso-compression	ratio end of shearing	iso-compression [kPa]	shearing [kPa]	strain [%]	strain [%]
	SPF1	0.593	0.513	0.454	300	500	8.9	20
	SPF2	0.407	0.368	0.382	300	500	1.7	32
	LBS1	0.544	0.476	0.460	1000	1700	7.2	20
	LBS2*	0.437	0.397	0.388	500	840	3.4	30
	LBS9*	0.551	0.490	0.466	500	950	5.5	23
	LBS10	0.436	0.319	0.281	5300	9700	8.6	17
	LBS14	0.435	0.405	0.422	100	190	0.9	21
	LBS16* 0.583	0.511	0.472	500	920	7.1	9
	LBS17* 0.457	0.409	0.386	500	1000	4.9	5
	LMS15 0.415	0.332	0.288	2400	5100	9.1	27

*selected samples for MIP and statistical analysis tests

Table 3

 3 Statistical parameters of SPF and comparable data of LBS and LMS.

	Test	Mean,	Median,	Standard	Skewness,	Kurtosis,	d60/d10 Diameter	Diameter
		μ [μm]	m	deviation,	γ	κ		at	at
			[µm]	σ [μm]				highest	smallest
								Peak	Peak
								[µm]	[µm]
	SPF-	0.76	0.71	5.77	0.63	4.42	10.36	1.60	0.58
	OED-LI
	SPF-	0.54	0.46	5.81	0.96	5.24	7.42	0.37	0.91
	OED-DI
	SPF-	0.57	0.47	6.18	0.97	4.90	7.40	0.42	1.02
	OED-
	LF
	SPF-	0.39	0.29	5.73	1.32	6.05	5.30	0.27	1.42
	OED-
	DF
	SPF1	1.15	0.52	11.89	0.77	2.62	8.30	-	-
	SPF2	0.68	0.46	7.85	1.05	4.37	7.40	-	-
	LBS2	1.85	1.89	4.49	0.78	5.86	6.01	-	-
	LBS9	1.84	1.98	3.78	0.88	6.35	6.75	-	-
	LBS16	1.84	1.79	4.12	0.81	5.99	5.35	-	-
	LBS17	1.79	1.63	4.12	1.22	6.59	4.79	-	-
	LMS-	2.71	1.76	4.78	1.88	4.21	3.07	-	-
	OED1
	LMS-	2.57	1.45	5.37	1.54	3.29	3.03	-	-
	OED2

ACKNOWLEDGEMENTS

The work described in this paper was fully supported by a grant from the Research Grants Council of the Hong Kong Special Administrative Region, China (project no. CityU 112813).

The authors are grateful to Prof. Pierre Delage for his valuable comments on MIP data. The technician Mr. Xavier Boulay and Baptiste Chabot are acknowledged for their technical support in conducting the MIP tests. Also, we would like to thank the reviewers for their constructive remarks which helped to improve the quality of the paper.

CONFLICT OF INTEREST STATEMENT

The authors certify that they have NO affiliations with or involvement in any organization or entity with any financial or non-financial interest (such as personal or professional relationships, affiliations, knowledge or beliefs) in the subject matter or materials discussed in this manuscript.

APPENDIX 1

Micro X-ray CT imaging

Visually, the micro X-ray CT scan images in Fig. A1 of dense and loose oedometer samples of LBS do not show significant differences in particle arrangement. The images also highlight a key problem for the use of CT to characterise the fabric of these soils. At the scale shown of 900μm across the image, there are insufficient larger particles and voids to be statistically reliable, and so much larger images would be needed. However, the sizes of the small particles and voids are already too small to be separated reliably by segmentation using existing methods within the resolution of the images, and this problem would be exacerbated for a larger image.

For this reason, MIP testing was preferred.