Calcium Environment in Silicate and Aluminosilicate Glasses Probed by 43 Ca MQMAS NMR Experiments and MD-GIPAW Calculations.

43 Ca MQMAS NMR spectra of three silica-based glasses in which Ca 2+ ions play different structural roles have been collected and processed in order to extract the underlying NMR parameter distributions. The NMR parameters have been interpreted with the help of molecular dynamics simulations and DFT-GIPAW calculations. This synergetic experimental-computational approach has allowed us to investigate the Ca environment, to estimate Ca coordination numbers from MD-derived models, and to push further the discussion about 43 Ca NMR sensitivity to the first and second coordination spheres: 43 Ca  iso and Ca-O distance can be successfully correlated as a function of Ca coordination number. Graphical abstract


Introduction
Ca is one of the most abundant elements in the Earth's crust and it is present in a wide variety of materials, either natural or manufactured, often playing a crucial role in the final physical-chemical properties of the system in which it is embedded.In silica-based glasses, depending on the compositions, Ca can play the role of network charge compensator or of network modifier.
For example, when added to silica glass, Ca ions disrupt two Si-O-Si bonds by forming two non-bridging oxygens (NBOs) leading to a less polymerized glass network and favouring the dissolution of the glass in water solution.On the contrary, in alumina-silicate and borosilicate glasses Ca ions can act as charge compensators of the excess negative charge of the [AlO4] -and [BO4] -units, and its consequently reduced mobility disfavours its dissolution in water [1].The situation is more involved when sodium ions are also present in the glass matrix since a competition between these ions can occur for occupying modifier of charge-compensator sites.[2] Since the role played by Ca affects glass macroscopic properties (viscosity, thermodynamic properties, conductivity, chemical durability, etc.) [3][4][5][6], the elucidation of its local environment, and therefore its structural role, is of fundamental importance for technological applications.
During the last decade, NMR spectroscopy has arisen as one of the leading edge techniques for the structural characterization of glasses.[7] Its high sensitivity to the local atomic environment, both in terms of chemical and topological features, stimulated the development of experimental protocols for solid state samples, like Magic Angle Spinning (MAS) [8,9] and two-dimensional (2D) multiple quantum (MQ) MAS experiments [10,11].
First 43 Ca NMR experiments started to appear from the late 90es, [12][13][14][15] and demonstrated the correlation between 43 Ca NMR parameters and Ca environment.However, 43 Ca solid state NMR is intrinsically difficult because 43 Ca is a quadrupolar nucleus (spin = 7/2) with a very low natural abundance (0.135%) and a small gyromagnetic ratio (γ).Thus it presents a very low NMR sensitivity with respect to other common, more abundant and high-γ nuclei, such as 23 Na [2,16] or 27 Al [17,18].To increase the signal to noise ratio, experiments at high magnetic fields with long acquisition times and/or isotope enrichment of samples are required; hence, the investigation of Ca environment is difficult and expensive [19].The combination of these factors results in a very limited number of NMR reports for the 43 Ca nucleus.So far, our knowledge about 43 Ca NMR experiments on silica-based glass samples is limited to two scientific works [19,20].
Shimoda et al. [20] reported very high-field 43 Ca MAS and 3Q/5Q/7QMAS NMR spectra of CaSiO3, CaMgSi2O6 and CaAl2Si2O8 with the aim to characterize Ca coordination number (CN).They distinguished three peaks in the 5Q and 7QMAS spectra and assigned them to three distinct Ca sites coordinated by 6, 7 and 8 oxygen atoms, with the latter two being the dominating CNs.From their analysis, calcium CN apparently increased from 7.1 in CaSiO3 to 7.8 in CaAl2Si2O8.However, the low sensitivity associated to the manipulation of high quantum transitions implies that they could have misinterpreted the 5Q and 7Q data, obtaining wrong distribution of coordination numbers.In fact, in a recent paper, some of us have demonstrated, through NMR computational spectroscopy [21], that a Ca coordination environment ranging from 5 to 8 CN well reproduced both the MAS and MQMAS spectra of CaSiO3 glass.[22] The sensitivity of 43 Ca NMR parameters to structural role played by Ca (charge compensator vs modifier) in silica and aluminosilicate glasses has been investigated by Angeli et al. [12] The 43 Ca MAS spectra of three glasses with composition 60SiO2•20Na2O•20CaO (CSN), 60SiO2•20Al2O3•10Na2O•10CaO (CASN) and 60SiO2•20Al2O3•20CaO (CAS) were collected at a magnetic field of 11.7 T and analysed by using the Gaussian Isotropic Model (GIM) for quadrupolar NMR parameters distribution.[23] Calcium acts as a charge compensator of [AlO4] -units in the CAS glass, as a modifier in the CSN glass and in both ways in the CASN glass.It was found that the transition from a charge-compensating role to network-modifying role results in an increase of both the isotropic chemical shift and the quadrupolar coupling constant of 43 Ca.For the CASN glass, the intermediate values of the NMR parameters suggest that calcium is probably in both types of environments.However, the CASN MAS spectrum could not be satisfactorily simulated using a combination of the spectra associated to CAS and CSN environment.
It was concluded that MQMAS experiments were needed to get more information.
In this work, we present the first 43 Ca 3QMAS NMR experiments of the samples investigated by Angeli et al. [12] The spectra have been analysed by using the Gaussian Isotropic Model (GIM) for fitting and an inversion method.On one hand the GIM method provides accurate information on the mean value and width of NMR parameter distributions, but it cannot reveal any asymmetry in the distribution if any, which shape is strongly influenced by the short-range environment of the nuclei.[2] On the other hand, with the inversion method the NMR parameter distributions are reconstructed numerically, without imposing any constraint on shapes and average values.The obtained experimental sets of NMR parameters have been compared to their theoretical counterparts obtained through plane waves Density Functional Theory (DFT) calculations with the Gauge Including Projector Augmented Wave (GIPAW) [24,25] approach on molecular dynamics (MD)-derived structural models, exploiting a recently introduced integrated computational approach.[21,[26][27][28][29] By coupling 43 Ca 3QMAS experiments, inversion procedure, and NMR-GIPAW calculations on sound MD-derived structural models of the glasses, we have been able to shed light on i) Ca short-range environment in terms of CN distribution and ii) the dependence of 43 Ca iso on Ca-O as a function of Ca CN improving our understanding of Ca environments in glasses.

Method
Experimental procedure CAS, CASN and CSN glass samples were synthetized by melting and quenching precursor oxides and carbonates.Details are provided in the work of Angeli et al. [12].In addition to original data, which consisted of MAS NMR spectra collected at a magnetic field of 11.75 T, here, 43 Ca MAS and 3QMAS NMR data were acquired on a Bruker 750MHz spectrometer, operating at a magnetic field of 17.6 T.
Samples were spun at 12.5 kHz and the chemical shift was measured respect to CaCl2(aq) 1M set at 0 ppm.For 3QMAS, a shifted-echo pulse sequence was used.[30] Typically 8192 signals were accumulated with a recycle delay of 1s for each t1 point.16 rotor-synchronized t1 increments (80 us) were collected for each 2D spectrum and processed using an in-house written software.
The idea underlying the inversion method is to reconstruct the NMR parameter distribution (  , ,   ) numerically, i.e. without using any analytical model.[17,31] But, there are some limitations because 43 Ca MQMAS spectra present smooth lines, i.e. without singularities arising from a well-defined or narrow ranged  values (sharp patterns), and are two-dimensional data.Hence, (  , ,   ) cannot be reconstructed numerically in three dimensions.As done before, [2,16,17] and as supported by MD simulations, an analytical model such as the marginal distribution of the GIM model for Q (Eq. 1) can be employed: By so doing, (  , ,   ) = (  ,   ) × () and only (  ,   ) needs to be reconstructed numerically, here with a Tikhonov regularization scheme [32] by imposing positivity and smoothness of (  , ,   ).[2] Sum of squared errors associated to this fitting procedure are actually very low (Table 1).

Computational details
For each composition, three independent structural models of 250 atoms have been generated by means of classical molecular dynamics (MD) using the DL_POLY package.[33] We used a force-field based on the core-shell model which has been successfully employed to simulate the structure of several silica-based glasses.[27,[34][35][36] For a full description of the forcefield functional forms and the parameters used in this work see ref. [27] and table S.1 of the ESI.
The glass compositions, numbers of atoms for each atomic species and densities are reported in table S.2 of the supplementary information (SI).The initial configurations were generated by placing randomly the atoms in a cubic cell which edge length was calculated using experimental densities [37] for CAS (14.98 Å) and CSN (14.98Å) glasses and the density estimated by the Priven empirical method [38] as encoded in Sciglass package [39] for the CASN glass (14.97Å).Each glass structure was generated according to the melt-quench approach in the NVT ensemble: the system was holded at 3200 K for 100 ps and cooled at a nominal rate of -5 K/ps to 300 K, then kept for 200 ps.[40] Integration of the equations of motion was run by the leap-frog algorithm every 0.2 fs; such a small value of the timestep is required to account for high frequency vibrations of the core-shell systems.[41] Once generated, final glass structures have been optimized by using plane-waves density functional theory (DFT) calculations with the CASTEP code [42].Constant pressure geometry optimization and a subsequent calculation of NMR parameters were performed at the Γ point [25], setting the convergence parameter for the Self Consistent Field energy as 5•10 -5 eV/atom, and the kinetic energy cutoff for the expansion of plane waves basis set as 700.0 eV.[22,43] The generalized gradient approximation (GGA) PBE functional [44] was employed and the corevalence interactions were described by ultrasoft pseudopotentials generated on the fly (see Gambuzzi et al. [27] for details).
Outputs of DFT-GIPAW calculations were processed by the fpNMR package, an in-house made code [22] devised to perform statistical analysis and to simulate NMR spectra (MAS, MQMAS, etc.).
Isotropic chemical shifts (δ iso calc ) were evaluated from the calculated isotropic magnetic shielding,    , using    = −(   −   ).For 43 Ca nucleus,   is 1134.1 ppm [22].Quadrupolar coupling constants, CQ, is calculated by employing the experimentally determined quadrupolar moment, Q, of -40.8 mB.[45] The CQ constants were also computed by using -44.4 mB [46] as suggested by Brugees et al. [47] in the case of Ca metal carboxylate compounds but we did not observe any improvement on the results.
Cutoff radii employed for the structural analysis (3.0 Å for Ca-O, 5.0 Å for Ca-Na and 5.0 Å for Ca-Al) were obtained from first minima of the radial pair distribution functions averaged over the last 201 configurations of the MD trajectories at 300 K.The coordination number of Calcium was calculated on the DFT-relaxed models according to the abovementioned cut-off radius.

Results
Figure 1 shows the high field 43 Ca experimental 3QMAS and MAS spectra of the CAS, CASN and CSN glasses.On one hand, shape MAS spectra (Figure 1.d) and isotropic projections, (Figure 1.e) are rather similar but widths are slightly different, suggesting similar but not identical ordering of Ca short-range environments in the three compositions.On the other hand, their positions are different, as observed in the previous study by Angeli et al. [12], denoting that, on average, a different role is played by Ca ions in the three glasses.The simulated spectra using the GIM and the inversion methods are also reported in Figure 2, from which it is evident that the second method better reproduces the tails of the experimental spectra.It is worth to note that a single line is observed on the MQMAS spectra for all compositions, especially for the CASN glass.
Figure 2 shows both experimental and theoretical NMR parameter distributions (CQ;iso) (panels g, h and j).The former are extracted with the inversion method from experimental MQMAS spectra, the latter results from MD-GIPAW calculations.Since the agreement of the two counterparts is very good for CAS and CASN and good for CSN, we can conclude that these models reproduce the vast majority of the chemical environments experienced by Ca in our samples.A further confirmation of the soundness of the MD-derived models is provided by the fact that theoretical MD-GIPAW 43 Ca 3QMAS spectra (simulated with Kernel Density Estimation formalism [22]) fit satisfactorily the experimental spectra, as shown in panels (k), (i) and (l).
The distributions of 43 Ca isotropic chemical shift ( 43 Ca iso) for the CAS, CASN and CSN glasses extracted by using both the GIM and the inversion fitting methods are displayed in Figure S1 of the ESI.The width of 43 Ca iso distributions are rather similar.In summary, 43 Ca CQ does not provide valuable information on the role played by Ca.This conclusion is also supported by the variation of CQ values with composition (Table 1).The computed CQ, are in agreement with experimental ones since they are within the experimental range of values determined by fitting 43  However, an in-depth analysis of the MD-derived structural features is justified by the good agreement between MD-GIPAW and experimental distributions of 43Ca isotropic chemical shifts.
This also allows us to investigate NMR/structural relationships.

Discussion
In the past, in analogy to the correlations found between 23 Na iso and <Na-O> in silicate crystals, [48,49] both experimental and computational works presented correlations between 43 Ca iso and <Ca-O> [14,15,[50][51][52] or Ca coordination number (CN) [51,52] in crystalline samples.However, in 2010 some of us [22] studied the CaSiO3 glass and demonstrated that 43 Ca iso barely correlates with Ca CN and <Ca-O> distances, and that only qualitative trends can be pointed out.Beside the lack of effective correlations, experimental uncertainties usually make NMR parameter distributions very close in shape and in width, making the associated structural information difficult to extract.
Nonetheless, mean values of experimental iso distributions extracted with the inversion fitting method (Figure S1  Amorphous computational models provide a wide variety of chemical environments experienced by 43 Ca, thus a higher number of points with respect to crystals.Under these conditions, the reliability of the aforementioned correlations can be verified.As in a previous work, [22] ambiguous trends between 43 Ca CQ and Ca first coordination shell (data not reported) are observed, even though 43 Ca CQ tends to slightly decreases with CN.
However, in line with the results obtained in a previous work for 23 Na, [2] the correlation between iso and <Ca-O> substantially improves when CN is expressly taken into account (Figure 3).Unluckily, the signal/noise ratio of the acquired 43 Ca 3QMAS spectra do not allow us to fit the experimental data and to obtain an estimation of the population of different Ca CN sites, as recently proposed Gambuzzi et al. [2] for Na in glasses.This impedes us to evaluate <Ca-O> by using a bivariate correlation.
However, MD simulations provide <Ca-O> and reveals that it becomes shorter moving from CAS (2.53 Å) to CSN (2.43 Å), showing an intermediate values in CASN glass (2.50 Å).This trend is reflected by an increase of the iso from CAS to CSN passing through the CASN glass.
These results are in good agreement with the ones reported by Angeli et al. [12], who estimated a <Ca-O> bond length of 2.58 Å, 2.47 Å and 2.36 Å for CAS, CASN and CSN glasses by exploiting an empirical correlation between 43 Ca iso and <Ca-O> elaborated on crystalline samples [15] and reported as a blue line in Figure 5.
In the latter, the intersection between the vertical black dashed lines, representing the values of 43  The quantification of the nature of the Ca role in the three glasses can be given by the ratio  = °  ⁄ , which ranges between 0 and 1.When X = 1, all O around Ca are BO, thus Ca plays the role of a pure network charge compensator; when X = 0 all O around Ca are NBO, thus Ca plays the role of a pure network modifier.showed that, when Na and Ca are both present, Na prefers to act as network charge compensator when Al is present in the network.This behaviour was also pointed out in previous computationalexperimental works [53,54].Nonetheless, extended Ca/Na intermixing around NBOs in Ca-Na silicate, aluminosilicate and borosilicate glasses was proved by both experimental [55] and computational [29] approaches.

Conclusions
In conclusion, the acquisition of 43 Ca MQMAS spectra allowed us to obtain, for the first time, a good set of NMR experimental data of the challenging 43 Ca nucleus in three silicate glass samples.Full iso and CQ distributions were provided by the inversion fitting procedure with small fitting errors.The set of theoretical NMR parameters, calculated on sound structural models of the glasses, provided us a valuable data set to find out that 43 Ca iso better correlates with <Ca-O> distances when CN is expressly taken into account, and to conclude that CN does not correlate with any NMR parameter.
Nonetheless, the study of Ca coordination in glasses has been possible thanks to the structural models: CN distribution is sharper when Ca acts as charge compensator and wider when it acts as network modifier.

Supporting Information
Tables about structural analysis of MD-derived and DFT-relaxed structural models.This material is available free of charge.Correlation between <Ca-O> and 43 Ca iso elaborated on crystalline samples (blue).<Ca-O> from MD-derived models (red) and 43 Ca iso obtained in this work from MQMAS spectra (dashed lines).For MD data, error bars represents width (standard deviation value) of the distribution.

Ca 3QMAS spectra presented in Figure 1 .
Both comp CQ and exp CQ variations are smaller than widths of their experimental and theoretical distribution (1.1-1.2MHz), thus making any assumption about the presence of trends not well founded.For CSN glass, a larger discrepancy is noted (despite a similar qualitative agreement from examination of figure2).This glass being richer in sodium cation known to be a more mobile species [REF] than Ca, the lower Cq experimental value could be attributed to thermal averaging effects [REF], as already recently noticed for 23Na [REF].Further investigation of such effects will necessitate finite temperature calculations (such as averaging over MD trajectories at room temperature), out of the scope of the present work.

Figure 4
Figure 4 reports the populations of the different X values in the three glasses.As expected, the

Figure 2 .
Figure 2.43 Ca 3QMAS experimental spectra (solid lines) for CAS, CASN and CSN glasses are compared to spectra simulated by using the GIM (panels a, b, c) and the inversion (panels d, e, f) fitting procedures.Panels g, h, j report the comparison of the experimental (dotted lines) and MD-GIPAW(full circles) (CQ;iso) distributions.Panels k, i, l report the comparison of experimental and MD-GIPAW (dashed lines)43 Ca 3QMAS spectra.

Figure 3 .
Figure 3. Theoretical 43 Ca iso versus <Ca-O> bond distances as a function of Ca CN for CAS (a), CASN (b) and CSN (c) glasses.Trend lines are presented in coherent colours for most populated CN sites.

Figure 4 .
Figure 4. Distribution of the ratio X = n°BOs/CN of Ca ions for the CAS, CASN and CSN glasses.

Figure 5 .
Figure 5. Correlation between <Ca-O> and43 Ca iso elaborated on crystalline samples (blue).<Ca-O> from MD-derived models (red) and43 Ca iso obtained in this work from MQMAS spectra (dashed lines).For MD data, error bars represents width (standard deviation value) of the distribution.
of the ESI) evidence that Ca experiences different short-range environments, which can be ascribed to different CN distributions.A confirmation derives from MD-derived models, where Ca CN is 5.7, 5.8 and 5.3 for CAS, CASN and CSN glasses, respectively (Table2).
Ca iso extracted in this work from MQMAS spectra, and the blue lines provides an estimation of the average <Ca-O> distances in the glasses (2.54 Å, 2.46 Å, 2.40 Å for CAS, CASN and CSN respectively) which are very close to the values provided by MD-derived models.The variation of <Ca-O> with composition has to be ascribed to the different roles played by Ca in the compositions under study: network modifier in CSN, charge compensator in CAS, both roles in CASN.As a network modifier, Ca depolymerizes the network and is coordinated mainly by NBOs, while when acting as network charge compensator Ca is allocated into polymerized network regions that need charge compensation for [AlO4/2] -, and is coordinated mainly by Bridging Oxygens (BOs).

Table 1 .
Average43Ca CQ and iso fitted from 3QMAS spectra ( exp ) and calculated on MD-derived models ( comput ), fitting errors (see SI for definition) and standard deviations on ( comput ) values are in parenthesis.