Ab Initio Thermodynamic and Thermophysical Properties of Sodium

Article pubs.acs.org/JPCA

Ab Initio Thermodynamic and Thermophysical Properties of Sodium Metasilicate, Na2SiO3, and Their Electron-Density and Electron-Pair-Density Counterparts Donato Belmonte,† Carlo Gatti,‡ Giulio Ottonello,*,† Pascal Richet,§ and Marino Vetuschi Zuccolini† †

DISTAV, Università di Genova, Corso Europa 26, 16132 Genova, Italy CNR-ISTM, c/o Dipartimento di Chimica, Università degli Studi di Milano, Via C. Golgi, 19, 20133 Milano, Italy § Institut de Physique du Globe de Paris, 1 rue Jussieu, 75005 Paris, France ‡

ABSTRACT: Thermodynamic and thermophysical properties of Na2SiO3 in the Cmc21 structural state are computed ab initio using the hybrid B3LYP density functional method. The static properties at the athermal limit are first evaluated through a symmetry-preserving relaxation procedure. The thermodynamic properties that depend on vibrational frequencies, viz., heat capacities, thermal expansion, thermal derivative of the bulk modulus, thermal correction to internal energy, enthalpy, and Gibbs free energy, are then computed in the framework of quasi-harmonic approximation. Acoustic branches are computed by solving the Christoffel determinant and are assumed to follow sine wave dispersion when traveling within the Brillouin zone. The procedure generates several thermo-physical properties of interest in materials science and geophysics (transverse and longitudinal wave velocities, shear modulus, Young modulus, Poisson ratio) all consistent with experimentally determined properties. A representative cluster is then abstracted from the cell and a detailed electron localization/ delocalization analysis is performed on it, in the ground state geometry, and on deformed states imposed by two peculiar mixed asymmetric stretching/bending modes affecting the silicate chain that, according to literature data, have anomalous mode Grüneisen parameters. A Bader analysis reveals an intriguing feature associated with these deformations: an increase in the covalence of the Si−O bond that strengthens the linkage opposing the weakening induced by thermal stress. Finally, on the same cluster, the Ramsey contributions to the JNM coupling are evaluated by the gauge-independent atomic orbital method. The calculated isotropic chemical shifts of both 23Na and 29Si are again in substantial agreement with observations. averaging of the biaxial chemical shift powder pattern in 29Si NMR spectra, points to rapid diffusion of Na+ near the melting point, inducing extensive librational motion of SiO4 tetrahedra.3 Through linear combination of atomic orbitals (LCAO) DFT/B3LYP ab initio procedures we calculated structure, elastic, and vibrational properties of crystalline Cmc21 Na2SiO3 in the aim of establishing a reference frame for the interpretation of the complex phenomena affecting the substance from 850 to 1362 K. The contribution of the acoustic modes was also evaluated along with the single-crystal elastic anisotropy. Clusters of the substance at variational state were then isolated to perform a detailed investigation on some molecular orbitals of particular relevance. A detailed analysis of the localization/delocalization indexes was performed to discover the points of weakness of the substance when subjected to thermal stress.

1. INTRODUCTION Sodium metasilicate (Na2SiO3) owes its practical importance to the fact that it is the crystalline form of the main component of water glass, of which millions of tons are produced annually as sources of silica or alkalis in chemical processes, as precursors for silicate synthesis or as binders. From a fundamental standpoint, this compound has the great interest of exhibiting a large premelting temperature interval of about 160 K over which various configurational changes take place in the solid state.1−3 These changes manifest themselves as premelting phenomena, which have been observed for a variety of silicates as marked excess heat capacities in temperature intervals ranging from a few tens to more than 200 deg below the congruent melting point.1 At the calorimetric onset of premelting, configurational changes set in in Na2SiO3, as observed by Raman and 29Si NMR spectroscopy where local rearrangements of the Si2O6 chains take place within a phase that still keeps its long-range order.2,3 The structure of Na2SiO3 was first determined by Grund and Pizy4 and later refined by McDonald and Cruickshank.5 More recently, a phase transition has been observed by both Raman spectroscopy and X-ray diffraction, which was interpreted as a transition from the Cmc21 space group to the Pmc21 space group.2 Moreover, Na site hopping, as deduced from the hightemperature Raman and 23Na NMR spectra and from partial © XXXX American Chemical Society

2. METHODS: STATIC CALCULATIONS AND EQUATION OF STATE (EOS) The primitive cell of Cmc21 sodium metasilicate (noncentrosymmetric) belongs to the orthorhombic pyramidal crystal class. Received: August 29, 2016 Revised: October 17, 2016 Published: October 19, 2016 A

DOI: 10.1021/acs.jpca.6b08676 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A It has four irreducible atoms in the asymmetric unit and 12 atoms in the unit cell. At the LCAO DFT/B3LYP level of theory,6,7 adopting the triple-ζ valence basis set with polarization of Peintinger and co-workers8 for Si and O and the 8-511G basis set of Dovesi and co-workers9 for Na, the structure rapidly converges to a stationary variational state with a0 = 10.5676 Å, b0 = 6.0977 Å, and c0 = 4.8523 Å. The optimized volume of the crystallographic unit cell is V0 = 312.673 Å3. Full geometry optimization (cell-edges and internal coordinates), was carried out with a symmetry-preserving relaxation procedure fully automatized in the CRYSTAL14 code.10 The agreement between the calculated and the experimental structural parameters is satisfactory (Table 1). The B3LYP volume is 2%

Table 2. Ab Initio Relative Axial and Volume Compression Data of Sodium Metasilicate (Na2SiO3), in the Pressure Range P = 0−3.5 GPaa

Table 1. Variational Structure (B3LYP) of Cmc21 Na2SiO3 at the Athermal Limit, Compared to Experimental Structural Refinement at Room Conditions

3. METHODS: VIBRATIONAL CALCULATIONS Optic Modes. The calculation of the phonon frequencies at Γ point was performed within the harmonic approximation diagonalizing the central zone (k = 0) mass-weighted Hessian matrix: Hij Wij(k ⃗=0) = Mi M j (1)

parameter

Na2SiO3

Na2SiO3

a0 b0 c0 a0/b0 c0/a0 Si−Obridging, Å Si−Ononbrdg, Å Si−O−Si (deg) Vcell, Å3 T, K ref

10.48 6.07 4.82 1.7265 0.4599 1.672 ± 0.002 1.592 ± 0.002 133.7 306.6 298.15 ref 5

10.5676 6.0977 4.8523 1.733 0.4592 1.689 1.619 133.042 312.673 0 this work

P (GPa)

a/a0

b/b0

c/c0

V/V0

0.0 2.2 5.4 9.1 13.5

1.0000 0.9915 0.9810 0.9707 0.9606

1.0000 0.9897 0.9775 0.9661 0.9561

1.0000 0.9897 0.9768 0.9628 0.9473

1.000000 0.971218 0.936646 0.902910 0.869986

a a/a0, b/b0, c/c0, and V/V0 are unit cell edges and volume ratios wih respect to reference values at the athermal limit (i.e., a0 = 10.568 Å, b0 = 6.098 Å, c0 = 4.852 Å, and V0 = 312.673 Å3).

where Mi and Mj are the atomic masses associated with the i,j-th Cartesian coordinates and Hij is the second derivative of the potential energy, evaluated at equilibrium (u = 0), with respect to the atomic displacements: ⎡ ∂ 2V ⎤ (x) ⎥ Hij = ⎢ u ∂ ∂ ⎢⎣ i uj ⎥⎦ 0

(2)

In eq 2, V(x) defines the potential energy surface (PES) of a periodic system with N nuclei (which is function of vector x of the 3N atomic coordinates), whereas ui = xi − x*i and uj = xj − x*j define the atomic displacements from the equilibrium positions (xi* and xj*). The optimized structure was taken as reference in the vibrational frequencies calculation. In the CRYSTAL code, first derivatives of the energy with respect to the atomic positions are calculated analytically for all uj coordinates,14 whereas second derivatives at u = 0 (where all first derivatives are zero) are calculated numerically.15 As already stated, there are 12 atoms in the Cmc21 orthorhombic unit cell of crystalline Na2SiO3 (i.e., Z = 2), so that there are 33 optic modes. The irreducible representations of the optic modes at Γ point are

higher than the experimental one, a well-known effect of exchange−correlation functionals based on generalized gradient approximation (GGA).11,12 Fitting the ab initio potential well with a third-order finite strain Birch−Murnaghan equation of state (BM3-EOS) (Figure 1) yields a static bulk modulus

Γoptic = 9A1 + 8A 2 + 9B1 + 7B2

(3)

All modes are Raman-active, whereas a total of 25 IR-active modes (9A1 + 9B1 + 7B2) are expected, along with three additional acoustic modes (1A1 + 1B1 + 1B2) corresponding to rigid translations of the lattice. In Table 3 we list the vibrational frequencies obtained in this study and their mode-γ analysis, which we carried out by repeating phonon calculations at various compressed states.16 The mode Grüneisen parameters (γi) define the volume dependence of the vibrational modes of the lattice in the framework of the quasi-harmonic approximation (QHA), being expressed as

Figure 1. Ab initio B3LYP potential well of Na2SiO3 (Cmc21). The adopted basis sets are TZVP for Si and O8 and 8-511G for Na9. The dashed line is an interpolation of the calculated energy points.

K0 = 70.81 GPa and a pressure derivative K′0 = 4.48. As far as we know, this is the first determination made so far and there are no other computational nor experimental data to compare with, except for some extrapolated data on polycrystalline phases along the join CaSiO3−Na2SiO3.13 Na−metasilicate displays anisotropic linear compressibilities with maximum compressibility along the c axis. The ratios a/a0, b/b0, and c/c0 predicted for Na2SiO3-Cmc21 in the pressure range P = 0−13.5 GPa are reported in Table 2.

∂ ln νi (4) ∂ ln V where νi is the wavenumber of the vibrational frequencies (expressed in cm−1). The ab initio IR and Raman spectra of Na−metasilicate at various T are shown in Figures 2 and 3. The infrared absorption intensities Πn were computed for each n-th mode by means of γi = −

B


Article

The Journal of Physical Chemistry A

in their intensities above ∼300 cm−1. As far as we know, these are the first determined IR spectra of the substance under investigation. In Figure 3a,b the ab initio Raman spectra at the athermal limit are compared with the experimental observations2 at room conditions. A close examination of Table 3 indicates that all experimentally observed Raman-active modes are reproduced by the calculation, though some of them do not appear explicitly in the spectra plotted in Figure 3 due to their very limited intensity. Like the IR modes, the Raman modes decrease progressively with increasing temperatures because all the relative mode Grüneisen parameters are positive (Table 3 and Figure 4). Acoustic Modes and Elastic Properties. At the longwavelength limit corresponding to the center of the Brillouin zone (k → 0), the first three solutions of the dispersion relation vanish.23 These solutions are related to the acoustic modes of the crystal and correspond to homogeneous translation of all the atoms in the structure along the three spatial directions by the effect of a polarized elastic wave. We obtained the phase velocity and polarization of the seismic wave along a given propagation direction by solving the Christoffel determinant:24

Table 3. Ab Initio B3LYP Vibrational Frequencies (νcalc, cm−1) and Relative Mode Grüneisen Parameters (γi) Calculated for Na2SiO3 Sodium Metasilicate (IR- and Raman-Active TO Modes), As Compared with the Experimental Data2 (νexp)a N

mode

νcalc (cm−1)

γi

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36

acoustic acoustic acoustic A2 A1 B2 A2 A1 A2 B1 B2 A2 A1 B1 B1 A2 B2 B2 A1 B2 B1 A2 A1 A2 B1 A1 B2 B1 A1 B1 A1 B1 B1 A1 A2 B2

87.5 103.2 167.8 142.2 169.4 179.7 184.1 189.9 192.5 192.8 209.6 223.4 227.2 229.2 239.0 246.0 247.3 259.7 278.9 301.9 304.3 387.2 392.1 409.4 415.1 497.3 528.6 544.7 575.4 700.3 866.7 876.1 948.6 954.8 1010.6 1012.2

1.5 1.5 1.5 0.4 1.6 2.2 0.8 3.6 3.0 4.0 2.3 2.2 3.8 1.8 2.4 3.4 1.0 4.0 1.9 1.7 1.1 0.7 1.0 0.4 0.8 0.4 0.2 0.3 0.5 0.5 0.6 0.5 0.5 0.5 0.6 0.6

νexp (cm−1)

Δcalc−exp (cm−1)

167

2

183

1

234

−5

|cijklnjnl − ρV 2δik| = 0 263

−3

302

0

394

−2

508

−11

549 589

−4 −14

868

−1

where V is the phase velocity of the seismic wave, n = n (nx, ny, nz) is the unit vector normal to the wave surface defining seismic propagation into space, cijkl is the elastic constant tensor, ρ is the density of the crystal, and δik is a Kronecker delta. The solutions are of two types: a quasi-longitudinal wave with polarization nearly parallel to the direction of propagation (VP), and two quasi-shear waves with polarization nearly perpendicular to ni (VS1 and VS2, the former propagating faster than the latter, i.e., VS1 > VS2). The components of the fourth-rank elastic constant tensor (cijkl) were computed using stress−strain relations as the second derivatives of the energy with respect to the strain components, that is cij =

965

−10

1015

−3

Acoustic frequencies are those calculated at the Brillouin zone boundary (kmax).

the mass-weighted effective mode Born charge vector,17,18 evaluated through a Berry phase approach.19,20 The isotropic Raman intensities Πn were computed for each mode through a fully analytical, coupled perturbed Hartree−Fock/Kohn−Sham scheme.21,22 The graphical representation of the infrared and Raman spectra S(ν) was obtained as a superposition of Lorentzian functions F, one for each mode:

∑ F(ν ;ν0,n ,Πn,δn) n

F(ν ;ν0, n ,Π n,δn) =

⎤ Πn ⎡ δn/2 ⎢ ⎥ π ⎢⎣ (ν − ν0, n)2 + δn 2/4 ⎥⎦

1 ∂ 2E V ∂εi ∂εj

0

(8)

In eq 8 we adopt a Voigt notation for the indexes (c11 → c1, c22 → c2, c33 → c3, c23 = c32 → c4, c13 = c31 → c5, and c12 = c21 → c6). We carried out the calculations by applying strains up to 0.02 in magnitude. The nine independent elastic stiffnesses (cij) of sodium metasilicate are listed in Table 4, along with the compliances (sij). Details about the stiffnesses−compliances conversion can be found elsewhere.25 We obtained the single-crystal seismic velocities (longitudinal and transverse) by solving the Christoffel determinant. However, these depend upon the relative orientation of the crystal with respect to the incident wave.26 Directionally averaged shear and longitudinal velocities may be then evaluated through the following scheme:27

a

S(ν ) =

(7)

u3 = VP,VRH ̅

(5)

u2 =

(6)

where δn is the damping factor of the n-th mode, which is related to the phonon lifetime. Being unable to compute this quantity, we used a constant value of 10 cm−1. S(ν) curves were evaluated in the range 0−1200 cm−1, in steps of 1 cm−1. Concerning the computed IR spectra (Figure 2a,b) all wavenumbers progressively decrease with increasing temperature without significant changes

u1 =

3

3

(9.1)

3 VS,VRH ·VS,MAX 3 ̅ 3 2(VS,VRH + VS,MAX 3) ̅

(9.2)

3 VS,VRH ·u2 3 ̅ 3 (2u 2 3 − VS,VRH ) ̅

(9.3)

where VS,MAX is the fastest shear-wave velocity in the single crystal, whereas V̅ P,VRH and V̅ S,VRH are the aggregate Voigt−Reuss−Hill C


Article


Figure 2. Computed IR spectra at selected temperatures: (a) low- to medium-frequency range; (b) high-frequency range.

longitudinal and shear seismic velocities28 (VRH), obtained from the aggregate elastic moduli (bulk modulus and shear modulus) as follows: KV =

1 × (c′ + 2c″) 9 −1

KR = (s′ + 2s″)

μR = 15 × (4s′ − 4s″ + 3s″′)−1

(10.1)

μ ̅VRH =

1 (μ + μ R ) 2 V

VS,VRH = ̅

μ ̅VRH /ρ

(10.2)

VP,VRH = ̅

K̅ VRH

1 = (KV + KR ) 2

(10.3)

μV =

1 × (c′ − c″ + 3c″′) 15

(11.1)

K̅ VRH +

(11.2)

(11.3) (12)

4 3

ρ

× μ ̅VRH (13)

In eqs 10.1 and 10.2 and 11.1 and 11.2, KV (and μV) and KR (and μR) are the Voigt and Reuss bounds for bulk modulus (and shear modulus), respectively; ρ is the density; c′ = c11 + c22 + c33 and D


Article


Figure 3. Computed Raman spectra (athermal limit): (a) low- to medium-frequency range; (b) high-frequency range. Experimental data of Richet and co-workers2 at 298.15 K (intensity and wavenumber) are superimposed for comparison (see also Table 3).

In this case we computed K̅ VRH = 70.74 GPa, compared to K0 = 70.81 GPa. Young’s modulus (E) and Poisson’s ratio (νP) are, respectively,

c‴ = c44 + c55 + c66 are obtained by summation of the diagonal terms of the stiffness tensor; c″ = c12 + c13 + c23 is the sum of the nondiagonal terms of the stiffness tensor; s′ = s11 + s22 + s33 and s‴ = s44 + s55 + s66 are obtained by summation of the diagonal terms of the compliance tensor; and s″ = s12 + s13 + s23 is the sum of the nondiagonal terms of the compliance tensor. As a proof of the internal consistency of the procedure, one must verify that the bulk modulus obtained by the analysis of the static potential well at the athermal limit is not much dissimilar from the aggregate bulk modulus obtained through the VRH averaging procedure (see Erba et al.29).

E=

νP = E

9K̅ VRHμ ̅VRH 3K̅ VRH + μ ̅VRH

(14)

3K̅ VRH − 2μ ̅VRH 2 × (3K̅ VRH + μ ̅VRH )

(15) DOI: 10.1021/acs.jpca.6b08676 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article


Figure 4. Ab initio B3LYP Raman spectra at various T conditions: (a) low- to medium-frequency range; (b) high-frequency range.

ν3 = 167.8 cm−1 (Table 3). All the relevant thermophysical properties calculated for crystalline Na2SiO3 (i.e., single-crystal, directionally averaged and aggregate elastic moduli and seismic velocities) are summarized in Table 4. Single-crystal azimuthal seismic anisotropy for longitudinal waves (AP) and shear waves (AS) may then be obtained as26

We obtained the maximum (angular) frequency of each acoustic branch (ωi, i = 1, 2, 3) from the directionally averaged seismic velocities (ui, i = 1, 2, 3 in eqs 9.1−9.3) by applying ⎛ 6π 2N ⎞1/3 0 ⎟ = uik max ωi = ui⎜ ZV ⎝ ⎠

(16)

where N0 is the Avogadro’s number, Z is the number of formula units in the unit cell, and V is the molar volume. The term kmax defines the Brillouin zone boundary (i.e., the radius af a sphere with the same volume of the Brillouin zone). Conversion of the obtained seismic velocities to linear frequencies yields the following results: ν1 = 87.5 cm−1; ν2 = 103.2 cm−1; F

⎛ VP,MAX − VP,MIN ⎞ ⎟⎟ × 100 AP = ⎜⎜ VP,VRH ̅ ⎝ ⎠

(17.1)

⎛ VS,MAX − VS,MIN ⎞ ⎟⎟ × 100 AS = ⎜⎜ VS,VRH ̅ ⎝ ⎠

(17.2) DOI: 10.1021/acs.jpca.6b08676 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article


frequency), and γi are the relative mode Grüneisen parameters as defined in eq 4. As for the product αKT of sodium metasilicate, it reaches a fairly constant value at high temperature, averaging 0.0047 GPa/K between 1000 and 2000 K. Following the guidelines developed in some of our previous works,12,31−35 we obtained

Table 4. Ab Initio B3LYP Elastic and Seismic Properties of Sodium Metasilicate (Na2SiO3)a c11 (GPa) c22 (GPa) c33 (GPa) c12 (GPa) c13 (GPa) c23 (GPa) c44 (GPa) c55 (GPa) c66 (GPa) K̅ VRH (GPa) K0 (GPa) μ̅VRH (GPa) V0 (cm3/mol) ρ0 (g/cm3) E (GPa) νP

s11 (GPa−1) s22 (GPa−1) s33 (GPa−1) s12 (GPa−1) s13 (GPa−1) s23 (GPa−1) s44 (GPa−1) s55 (GPa−1) s66 (GPa−1) V̅ S,VRH (Km/s) V̅ P,VRH (Km/s) VS,MIN (Km/s) VS,MAX (Km/s) u1 (Km/s) u2 (Km/s) u3 (Km/s)

136.032 112.798 130.928 50.912 35.334 42.654 59.282 41.086 20.666 70.74 70.81 38.44 47.073 2.593 97.59 0.2701

9.0253 × 10−3 11.5434 × 10−3 8.8880 × 10−3 −3.5955 × 10−3 −1.2644 × 10−3 −2.7903 × 10−3 16.8687 × 10−3 24.3389 × 10−3 48.3884 × 10−3 3.85 6.86 2.83 4.78 3.58 4.22 6.86

∂KT ∂T

of crystalline Na2SiO3 is reproduced by the following polynomial function: αT = α0T + α1 + α2T −1 + α3T −2 + α4T −3

3R ⎛⎜ 2 ⎞⎟ Z ⎝π ⎠

CV =

c11, c22, c33, c12, c13, c23, c44, c55, c66 = independent components of the elastic stiffness tensor; s11, s22, s33, s12, s13, s23, s44, s55, s66 = independent components of the elastic compliance tensor; K̅ VRH = Voigt−Reuss−Hill average of the aggregate bulk modulus; K0 = single-crystal bulk modulus (static value); μ̅VRH = Voigt−Reuss−Hill average of the shear modulus; V0 and ρ0 = static molar volume and density at the athermal limit (i.e., T = 0 K and P = 0 GPa), respectively; E = Young’s modulus; νP = Poisson’s ratio; V̅ P,VRH, V̅ S,VRH = Voigt−Reuss−Hill average of the aggregate shear and longitudinal seismic velocities, respectively; VS,MIN, VS,MAX = minimum and maximum single-crystal shear-wave velocities, respectively; u1, u2 = directionally averaged shear-wave velocities; u3 = directionally averaged longitudinal-wave velocities.

+

3n i

⎛

i=4

⎞2 ⎟ e Xi − 1 ⎠

3

3n i

⎛

4

(Xi 2 − X2)1/2 (e X − 1)2

0

i=1

∑ e X ⎝⎜

[arcsin(X/Xi)]2 X2 e X dX

Xi

∑∫

⎞2 ⎟ e Xi − 1 ⎠ Xi

(20)

CP = CV + TαT 2KP , T VP , T

(21)

we evaluated CP (Figure 6) and accounted for its variations between 298.15 and 2000 K by the polynomial expansion commonly adopted in metallurgy:38 CP = a + b × T + c × T −2 + d × T −1/2 + e × T −3 + f × T2 + g × T3

(22)

where a = 72.351, b = 48.796 × 10−3, c = 9.1136 × 106, d = 545.24, e = −29.913 × 108, f = 27.157 × 10−6, and g = −11.304 × 10−9. The calculated heat capacities show a very good agreement with the calorimetric results of Kelley39 at low T and with those of Richet and co-workers1 up to the limit of the premelting perturbative effect (see Discussion) (Figure 6). In the absence of configurational disorder, the entropy of one mole of a crystalline substance reduces to the following expression: 3

S=

4. THERMOCHEMISTRY Thermoelastic properties can be obtained from the quasiharmonic expression of the αKT product:

∑ γie X ⎜⎝

R Z

3

where n is the number of atoms in the unit cell and Z is the number of formula units per unit cell. The first term on the righthand side of eq 20 is the acoustic contribution at kmax represented in terms of a sine-wave dispersion relation.27,36,37 The second term is the contribution of all remaining (optic) modes. The isochoric heat capacity computed at the B3LYP level attains, as expected, the Dulong−Petit limit at high temperature. Using the standard relation between CP and CV, i.e.:

where VMAX and VMIN are the fastest and slowest longitudinal and shear velocities, respectively, and V̅ P,VRH and V̅ S,VRH are isotropic seismic velocities averaged over all propagation directions. The obtained results are presented in Figure 5 by means of inverse pole figures plotted on contoured stereograms using the 30 CAREWARE package and show that crystalline Na2SiO3 has relatively high AP (∼15%) and extremely high AS (∼51%) maximum anisotropies, hence a strongly anisotropic single-crystal elasticity. This is also confirmed by the elastic tensor calculation, whereby c11, c22, c33 values are much higher (85 GPa on average) than those of the other stiffnesses (Table 4). Concerning seismic velocities, the [011] and [110] crystallographic directions are those of fastest and slowest propagation for longitudinal waves (VP), respectively; c and a crystallographic axes are those of fastest and slowest propagation for shear waves (VS) (Table 4 and Figure 5). A notable fact is that the shear-wave seismic anisotropy displays low to moderate values (from ∼1% to 20%) throughout the volume of the crystal, except along the b axis, where it quickly rises to a huge maximum of ∼51%.

R ZV

(19)

The numerical coefficients obtained by regression of the ab initio data are α0 = 4.178 × 10−9; α1 = 2.1125 × 10−5; α2 = 2.3383 × 10−2; α3 = −6.5979; α4 = 550.379. In the harmonic approximation the isochoric heat capacity of one mole of substance is

a

αKT =

( barK ) (≈−0.017 GPa/K). The thermal expansion

= −168

3 3R ⎛⎜ 2 ⎞⎟ ∑ Z ⎝ π ⎠ i=1

−

Xi

+ (18)

+

where R is the gas constant, V is the molar volume, Z is the number of formula units in the unit cell, n is the number of atoms in the unit cell, Xi is the undimensionalized frequency of the i-th vibrational mode (Xi = hυi/kT, where h and k represent Planck and Boltzmann constants, respectively, and υi is the radial

3R ⎛⎜ 2 ⎞⎟ Z ⎝π ⎠ R Z

∫0

3n

⎛

∑ ⎝⎜

i=4 T

3

∫0

(Xi 2 − X2)1/2 (e X − 1)

3

∑∫ i=1

[arcsin(X /Xi)]2 X dX

Xi

Xi

0

[arcsin(X /Xi)]2 2

2 1/2

(X i − X )

ln(1 − e−X ) dX

⎞ R − ln(1 − e−Xi)⎟ + ln(Q e) ⎠ Z e −1 Xi

Xi

VTKT αT 2 dT

(23)

The first two terms on the right in eq 23 are sine-wave dispersion contributions (acoustic modes) of the Kieffer model,27,36,37 the third is the contribution of the remaining 3n − 3 (optic) modes, G


Article


Figure 5. Inverse pole figures of longitudinal seismic velocities (VP) and quasi-shear-wave velocities (VS1 and VS2) calculated for Cmc21 Na2SiO3 using the elastic tensors obtained in this study. Single-crystal seismic anisotropies (AP and AS, in %) and the difference between maximum and minimum shearwave anisotropy (AVS, in %) are also reported at the bottom of each stereogram. X1, X2, and X3 are the reference axes, with X1 along [100], X2 along [010], and X3 along [001] crystallographic directions. Contour lines are in (K m)/s.

Figure 6. Ab initio B3LYP isobaric heat capacity of Cmc21 Na2SiO3 (solid line) compared to experimental values1,39,42 and tabulated data.41

However, because α has been derived “quasi-harmonically”, the estimate of the last integral in eq 23 is a quasi-harmonic estimate as well. The bulk entropy, evaluated according to eq 23, attains at 298.15 K a value of 111.853 J/(mol × K), in reasonable agreement with the thermochemical tables.40,41 The overall

the fourth term is the electronic contribution arising from spin multiplicity (here identically equal to zero) and the last is the anharmonic contribution to the entropy of the substance. The anharmonic term arises formally from the differential increment of isobaric dT and isochoric heat capacities with T (i.e., dSanh = (CP − CV ) T ). H


Article


Figure 7. Ab initio B3LYP entropy of Cmc21 Na2SiO3 compared to thermodynamic tables.40,41

agreement over the entire T range here investigated is also satisfactory (Figure 7). In the harmonic approximation the internal energy of a solid is given by U = E B3LYP,crystal + +

RT Z

3n

3RT ⎛⎜ 2 ⎞⎟ Z ⎝π ⎠

⎛

∑ Xi⎜⎝ 1

2

i=4

+

3 3

∑∫ i=1

0

Xi

Table 5. Thermochemical Reference Data for the Various Atoms of Interest in This Studya

(24)

∑ niHf,A ,0 − ∑ niΔHelement ,0 → 298.15 i

i=1

ref

445.668 246.790 107.550.

3.2180 4.3415 6.4475

−289.390650 −75.095257 −162.224203

35 35 this work

H0f,0 = enthalpy of formation of the gaseous atom from the stable element at T = 0 K, P = 1 bar (NIST-JANAF Tables41). ΔHelement,i,0→298.15 is the enthalpy difference between T = 298.15 K and T = 0 K for the monatomic element (NIST-JANAF Tables41). EAi are the (optimized) electronic energies of the gaseous atoms. n/Z

D0 =

∑ niE A i=1

i

− (E B3LYP,crystal + EZPE,crystal )

(26)

This energy term corresponds to the sum of the electronic energies of the gaseous atoms at 0 K minus the electronic + zero point energy of the crystal. Concerning the electronic energies of the isolated atoms, it must be stressed that the basis sets adopted for atoms-in-crystals are not adequate to explore outer regions appropriate for the gaseous state. For example, the 8-511G basis set adopted for Na returns at the DFT/B3LYP theory level an electronic energy of E0 = −162.023898 hartree, which is intermediate between the B3LYP/STO-3G low limit of −160.082859 hartree and the B3LYP/Aug-cc-pVQZ value of −162.298993 hartree. Adopting optimized electronic energies for Si and O,35 the electronic energy of Na returning the experimental value of H0f is E0 = −162.224203 hartree (Table 5). This computation gives the correct dissociation energy of the substance (Table 6) if one assumes that the basis set superposition error (BSSE) is negligible. The thermodynamic properties and thermophysical parameters of crystalline Na2SiO3 as obtained in this study are summarized in Table 6. A detailed thermodynamic tabulation in the temperature range from 298.15 to 2000 K is reported in Table 7.

n/Z i

EAi (hartree)

Si O Na a

0 = (ΔU298.15 − 0 + P ΔV298.15 − 0) − D0 H̅ f,298.15

i=1

ΔHelement,i,0→298.15 (kJ/mol)

(Xi 2 − X2)1/2 (e X − 1)

where the second and third terms on the right constitute the thermal correction to U (ΔUT−0, including the zero-point energy). Knowing U and S, the absolute enthalpy of a substance (H) and the absolute values of its Helmholtz and Gibbs free energies (F and G, respectively) could be obtained by the usual thermodynamic relations (H = U + PV; F = U − TS; G = H − TS). These absolute values should not be confused with the enthalpy and free energy of formation of a substance from elements or oxides that are given in the tabulations of thermodynamic properties. Calculation of the ab initio enthalpy of formation of a substance from the elements requires the assessment of a thermochemical cycle of the type:

+

H0f,Ai,0 (kJ/mol)

[arcsin(X /Xi)]2 X dX

1 ⎞ ⎟ Xi e − 1⎠

n/Z

species

(25)

where the term within brackets is the thermal correction to the enthalpy at 298.15, Hf,Ai,0 is the enthalpy of formation of the i-th gaseous atom from the stable element at T = 0 K, P = 1 bar (with the summation extended to the n/Z atoms in the molecule), ΔHelementi,0→298.15 is the enthalpy difference between T = 298.15 and T = 0 for the monatomic elements (Table 5), and D0 is the zero-point dissociation energy of the gaseous molecule into gaseous atoms at 0 K: I


Article


tensor between sodium and oxygen in Cmc21 cluster [Na8Si4O13]2− are dominated by the interaction between the magnetic moments associated with the nuclear spins (i.e., the Fermi contact, FC). The JNM coupling tensor between sodium and oxygen is clearly affected by the interionic distance and decidedly more intense when nonbridging oxygen centers (nb in Table 8) are involved. The fact that all factors are reasonably similar (besides those ones involving terminal parts of the chain; b′ in Table 8) stresses again that the cluster is a good proxy of the crystalline structure. The calculated isotropic chemical shifts of 23Na with respect to the isotropic shielding of a one-molal NaCl solution computed at the same theory level through IEFPCM45 [σref = 622.8 ppm; δiso = −(σ − σref) = 33−51] is slightly higher than what observed in the crystalline substance3 (δiso ≈ 25) although substantially comparable when one considers that the experimental 23Na static spectrum is “broadened to the point of being featureless”.3 The calculated isotropic shift of 29Si with respect to the isotropic shielding of tetramethylsilane [Si(CH3)4] calculated at the same theory level [σref = 347.2 ppm; δiso = −(σ − σref) = −51 to − 75] is again in substantial agreement with experimental observations3 (δiso = −79.6). For the isobaric heat capacity, our first-principles calculations deviate from the calorimetric values only above ∼1000 K, where the Dulong−Petit limit is almost completely reached. Because an increase of the bulk modulus with T is unphysical, the anomalous increase in anharmonicity (i.e., TVα2K) causing the observed deviation can only be eventually ascribed to an anomalous thermal expansion (within the homogeneity range of the substance) not revealed by the γ-mode analysis. The fact that this anomalous behavior is precursor to melting suggests that melting itself is caused by the very same phenomenon when nuclear centers are too far to result in a sufficient electronic delocalization between adjacent basins. Further bond-strength weakening may arise from the loss of electrostatic (stabilizing) effect on outer bonding molecular orbitals (MO) of the [SinO3n+1]−2n−2 chain when cationic centers move away from the chain itself. In fact, as already discussed elsewhere,46 these effects, which may be precisely evaluated for inner orbitals in terms of Kohn−Sham core-level (O 1s) energy shifts47−49 are also observed on the low-lying σ-bonding and antibonding molecular orbitals (MOs) of the oxygen framework in zeolitic rings. Clear relationships with the core-level shifts are moreover found when the effects of counterions are evaluated through electron localization/delocalization indices (i.e., Bader analysis, see later on), which are defined in terms of the whole wave function and not just of the individual orbitals. Figure 10 shows the σ-bonding orbital of the [Na8Si4O13]2− cluster (MO N. 50, out of 273; 125 occupied). The electron cloud (drawn at an isodensity limit of 0.03 e/Å3) encompasses essentially all the atoms of the chain besides the terminal ones. By analogy with what was already seen for zeolitic rings,46 we interpret this MO as a Pauling double bond, contributing substantially to the structural stability of the chain. The eigenvalue of MO is −1.087 hartree, quite similar to the energy levels observed in similar MOs of zeolitic rings.46 When Na is replaced by H while keeping the symmetry and structure of the cluster constant, the electron cloud keeps the same shape, whereas its eigenvalue rises to −1.061 hartree. To further investigate the role of deformation of structural stability, we submitted the cluster to an electron localization/ delocalization analysis. Electron localization/delocalization is related to the two-electron, or pair, density, which contains all information about the correlated motion of a pair of electrons.50,51 Several descriptors, measuring the extent to which electrons are

Table 6. Ab Initio B3LYP Thermodynamic Properties of Sodium Metasilicate (Na2SiO3)a V0298 (cm3/mol) S0298.15 [J/(mol × K)] a b × 103 c × 10−5 d e × 10−8 f × 106 g × 109

47.326 111.85 72.351 48.796 91.136 545.24 −29.913 27.157 −11.304

D0 (kJ/mol) K0 (GPa) K′0 (dK/dT)P (bar/K) α0 × 109 α1 × 105 α2 × 102 α3 α4

2951.87 70.81 4.48 −168.0 4.178 2.1125 2.3383 −6.5974 550.379

a 0 V 298 = molar volume at standard state (T = 298.15 K and P = 1 bar); S0298.15 = standard state entropy; a, b, c, d, e, f, g = numerical coefficients of the polynomial function for isobaric heat capacity (CP), fitted in the T range 298.15−2000 K (eq 22); D0= dissociation energy at the athermal limit (see text); K0 = static bulk modulus; K′0 = pressure derivative of bulk modulus; (dK/dT)P = isobaric temperature derivative of bulk modulus; α0, α1, α2, α3, α4 = numerical coefficients of the polynomial function for thermal expansion (eq 19).

5. DISCUSSION When comparing the ab initio isobaric heat capacity of the substance (Tables 6 and 7) with the calorimetric measurements,42 we see that the premelting zone (from ∼1030 K to the melting point 1362 K) results in a limited increase in enthalpy and entropy (Figure 8a,b respectively). The second-order contribution to melting is, however, almost negligible in terms of Gibbs free energy (i.e., ∼0.6 kJ/mol at Tf = 1362 K). There is, in this respect, an interesting analogy to λ contributions in solid state phase transitions. Experiments and first-principles calculations are also in substantial agreement regarding the Raman spectra of the low-temperature Cmc21 phase. The mode Grüneisen parameters are all intrinsically positive (Table 4); hence, Raman modes are expected to decrease with T, as experimentally observed.2 However, the Raman-active mode ν12 (two-components ν12 and ν12a) ehibits a nonlinear shift with temperature above ∼900 K and bands ν12 and ν17 (assigned respectively to Si−O−Si bending + stretching and stretching motions at ν12 = 588 cm−1, ν17 = 965 cm−1) have a sudden increase in the area ratio with T about 100 K below the melting point of the substance.2 Mode ν12 corresponds to our ν29 A1 motion (575.4 cm−1; γi = 0.5) and mode ν12 to the ν34 A1 motion (954.8 cm−1; γi = 0.5). In the aim of investigating in detail these two peculiar vibrations, we abstracted from the variational structure a cluster corresponding to an elongated packing along the c axis. The cluster contains a tetrameric chain surrounded by Na atoms in the positions dictated by the Cmc21 symmetry (Figure 9). The spectral analysis carried out with the GAUSSIAN package43 at the same level of theory adopted for the crystal assigns to the band ν54 of the cluster, which corresponds to the ν29 mode of the crystal, a frequency of 597 cm−1. The mode is effectively a mixed Si−O bending + Si−O−Si stretching motion affecting mainly the tetrameric chain (see the displacement vectors in Figure 9a). The Na atoms are immobile and completely unaffected by the thermal stress. We show in Figure 9b the displacement vectors of the ν64 mode of the cluster (954.2 cm−1), which corresponds to the ν34 of the crystal (954.8 cm−1). This vibrational mode is essentially described by a stretching motion mostly affecting the bridging oxygens of the chain.2 Again, the Na atoms are unaffected by thermal vibration. As evaluated by the gauge-independent atomic orbital (GIAO) method44 (Table 8), the Ramsey contributions to the JNM coupling J


Article

The Journal of Physical Chemistry A Table 7. Ab Initio B3LYP Thermodynamic and Thermophysical Functions of Sodium Metasilicate (Na2SiO3)a T (K)

KT,P (GPa)

α × 105 (K−1)

Pthermal (Gpa)

CV [J/(mol × K)]

CP [J/(mol × K)]

S [J/(mol × K)]

Uvib (kJ/mol)

Hvib (kJ/mol)

Gvib (kJ/mol)

γthermal

298.15 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000

68.84 68.82 67.63 66.32 64.90 63.40 61.80 60.12 58.33 56.43 54.41 52.26 49.97 47.50 44.85 41.97 38.84 35.41 31.60

5.856 5.868 6.366 6.679 6.939 7.192 7.454 7.732 8.024 8.330 8.648 8.977 9.315 9.661 10.014 10.373 10.738 11.107 11.480

0.732 0.739 1.157 1.594 2.042 2.497 2.956 3.419 3.884 4.350 4.818 5.287 5.757 6.227 6.698 7.170 7.641 8.113 8.586

109.234 109.534 122.191 130.045 135.112 138.511 140.878 142.581 143.842 144.799 145.542 146.128 146.599 146.983 147.299 147.563 147.785 147.974 148.136

112.567 112.900 127.396 137.141 144.190 149.701 154.343 158.505 162.428 166.265 170.122 174.068 178.156 182.429 186.917 191.651 196.653 201.948 207.557

111.853 112.549 147.107 176.572 202.135 224.682 244.863 263.161 279.937 295.467 309.967 323.604 336.514 348.804 360.562 371.862 382.760 393.310 403.545

61.8 62.0 73.6 86.3 99.5 113.2 127.2 141.4 155.7 170.2 184.7 199.3 213.9 228.6 243.3 258.0 272.8 287.6 302.4

61.8 62.0 73.6 86.3 99.6 113.2 127.2 141.4 155.7 170.2 184.7 199.3 213.9 228.6 243.3 258.0 272.8 287.6 302.4

28.9 28.7 15.9 0.0 −18.5 −39.3 −62.0 −86.6 −112.7 −140.2 −169.0 −199.1 −230.3 −262.5 −295.8 −329.9 −364.9 −400.7 −437.4

1.747 1.746 1.673 1.634 1.614 1.605 1.603 1.605 1.610 1.618 1.627 1.638 1.651 1.664 1.679 1.694 1.711 1.729 1.746

KT = isothermal bulk modulus; α = thermal expansion coefficient; Pthermal = thermal pressure; CV, CP = isochoric and isobaric heat capacity; S = entropy; Uvib, Hvib, and Gvib = vibrational contributions to internal energy, enthalpy, and Gibbs free energy, respectively. a

Tf

ΔCP dT ; Figure 8. Excess enthalpy (a) and excess entropy (b) terms obtained by integrating the heat capacity difference Hexcess = ∫ 298.15

Sexcess =

T

f ∫298.15

ΔCP T

dT between the computational CP obtained in this study (Table 6) and the calorimetric results of Richet and co-workers.1

Figure 9. [Na8Si4O13]2− cluster, obtained from 1 × 1 × 3 packing of the Cmc21 variational structure of the crystal. The displacement vectors affecting Si and O atoms in the chain are superimposed: left side, ν54 = 597.1 cm−1; right side, ν64 = 954.2 cm−1.

shared between two or more atoms (or other chemical “entities”) have been introduced over the years.52 They are usually defined in terms of the exchange−correlation density (XCD), which

represents the deviation, caused by the Coulomb and Fermi correlation of electron motions, between the true pair density of a system and that given by the purely classical description of a K


Article


Table 8. Ramsey Contributions to the JNM Coupling Tensor between Sodium and Oxygen and between Silicon and Oxygen in the [Na8Si4O13]2− Cluster Evaluated with the GIAO Procedure at the DFT/B3LYP Theory Level and with the Same Basis Set Adopted for the Cmc21 Na2SiO3 Crystala couple

coordination

dM−O (Å)

FC

SD

PSO

DSO

JNM

Na24O10(nb) Na25O6(nb) Na22O5(nb) Na23O2(nb) Na19O12(b) Na21O9(b) Na18O8(b) Na20O7(nb′) Si15O10(nb) Si15O11(nb) Si15O12(b) Si15O13(b′) Si14O12(b) Si14O6(nb) Si14O9(b) Si14O4(nb) Si17O9(b) Si17O5(nb) Si17O8(b) Si17O7(nb′) Si16O8(b) Si16O2(nb) Si16O3(b′) Si16O1(nb)

VI VI VI VI ? ? ? ? IV IV IV IV IV IV IV IV IV IV IV IV IV IV IV IV

2.29 2.29 2.29 2.29 2.42 2.42 2.42 2.37 1.62 1.62 1.69 1.70 1.70 1.62 1.69 1.62 1.70 1.62 1.69 1.62 1.70 1.62 1.69 1.62

−20.9 −37.7 −33.6 −32.3 −6.5 −4.1 −3.7 −15.7 19.4 24.8 16.8 25.6 12.6 16.5 9.7 20.6 7.6 15.5 7.9 19.4 12.2 16.8 22.6 26.0

−0.5 −0.1 −0.1 0.0 −0.1 −0.1 −0.1 −0.1 0.0 −0.3 −0.4 −1.4 −0.5 −0.4 −0.4 −0.4 −0.4 −0.4 −0.4 −0.4 −0.5 −0.3 −0.4 −0.4

0.5 1.0 0.7 0.9 1.0 0.8 0.8 1.3 7.0 9.5 4.2 10.5 5.6 8.6 5.4 9.3 4.9 8.5 5.7 9.0 4.1 7.9 8.0 9.9

−0.1 −0.1 −0.1 −0.1 −0.1 −0.1 −0.1 −0.1 −0.1 −0.1 −0.1 −0.1 −0.1 −0.1 −0.2 −0.1 −0.2 −0.1 −0.1 −0.1 −0.1 −0.1 −0.1 −0.1

−21.0 −36.9 −33.0 −31.4 −5.7 −3.5 −3.1 −14.2 26.1 33.9 20.5 34.5 17.6 24.5 14.5 29.4 12.0 23.5 13.1 28.2 15.7 24.3 30.1 35.5

a dM‑O = metal to first neighbour distance; FC = Fermi contact; SD = spin dipolar interaction; PSO = paramagnetic spin−orbit term; DSO = diamagnetic spin−orbit term. Distances are in Å; all other values are expressed in Hz. See Figure 9 for atom labeling.

Figure 10. σ-bonding inner orbital of the [Na8Si4O13]2− (left side) and [H8Si4O13]2− (right side) clusters resolved at the B3LYP level of theory with Cmc21 symmetry and constant atomic coordinates. The adopted basis set is the triple-zeta valence basis set with polarization of Peintinger and co-workers8 for Si and O, the 8-511G basis set of Dovesi et al.9 for Na, and the 5-11G* set of Dovesi and co-workers58 for H. Electron clouds are drawn for an isodensity limit of 0.03 e/Å3. The eigenvalue for [Na8Si4O13]2− is −1.087 hartree, and that of [H8Si4O13]2− is −1.061 hartree.

(exchanged) between Ω and Ω′ (with Ω ≠ Ω′). In practice, N(Ω) is the sum of the number of electrons that, on a statistical basis, always stay in Ω and of half the number of those that are exchanged between Ω and all the remaining basins. The λ and δ indices lend themselves to insightful statistical interpretations of electron localization/delocalization effects because they are intimately related to the variances and covariances of the electron populations of their associated domains.54 For our Kohn−Sham DFT wave functions, we calculated the λ(Ω) and δ(Ω,Ω′) indices with the customarily adopted, yet approximate, formula25 that makes use of a Hartree−Fock-like XCD matrix. The indices determined in this way are known to be often very close to those obtained with the Hartree−Fock (HF)

product of independent electron densities. In the approach of Fradera and co-workers,53 the XCD is doubly integrated over either the same QTAIM (quantum theory of atoms in molecules) atomic domain, Ω, or over two different domains, Ω and Ω′, to yield the so-called localization [λ(Ω,Ω); hereinafter λ(Ω), for the sake of simplicity] or delocalization δ(Ω,Ω′) indices, respectively. The number of electrons N(Ω) on an atom Ω is then given by N(Ω) = λ(Ω) +

1 2

∑ Ω′≠Ω

δ(Ω′, Ω) (27)

where λ(Ω) is the number of electrons that are fully localized within Ω and δ(Ω,Ω′) the number of those that are delocalized L


Article

The Journal of Physical Chemistry A theory,54,55 though such similarity does not imply that an electron-pair density matrix may be obtained at the density functional level of theory with the formalism valid for HF wave functions. Table 9 lists the interionic distances of the [Na8Si4O13]2− cluster in the Cmc21 ground state geometry (GST) of Na2SiO3

Table 10. Electron Populations (Bader’s Net Atomic Charges) of Atoms Ωa

Table 9. Electron Localization/Delocalization Analysis of the [Na8Si4O13]2− System at the Ground Structural State (GS) and Subject to ν54 and ν64 Maximum Deformationa X−Y couple Si15−O13

Si15−O12

Si14−O12

Si15−O11

Si15−O10

O11−Na19

O10−Na24

Si15−Si14

R (Å) ρbcp, au δ(X,Y) R (Å) ρbcp, au δ(X,Y) R (Å) ρbcp, au δ(X,Y) R (Å) ρbcp, au δ(X,Y) R (Å) ρbcp, au δ(X,Y) R (Å) ρbcp, au δ(X,Y) R (Å) ρbcp, au δ(X,Y) R (Å) ρbcp, au δ(X,Y)

GS

ν54

ν64

1.6988 0.1286 0.4171 1.6892 0.1173 0.3090 1.6988 0.1220 0.3382 1.6188 0.1484 0.4718 1.6188 0.1446 0.4302 3.8511 no bcp 0.0012 2.2875 0.0293 0.2224 3.1075 no bcp 0.0030

1.8442 0.1016 0.5394 2.0132 0.0624 0.2345 1.9885 0.0800 0.4014 1.7063 0.1265 0.5993 1.7161 0.1202 0.5204 3.8437 no bcp 0.0006 2.1376 0.0374 0.2225 3.5930 no bcp 0.0021

1.7358 0.1198 0.3682 1.1994 0.3907 0.5802 2.0182 0.0698 0.3126 1.7916 0.1071 0.4050 1.7916 0.1027 0.3591 3.89478 no bcp 0.0022 2.23075 0.0339 0.2864 2.87209 no bcp 0.0048

Ω

GS

ν54

ν64

Si14 Si15 O10 O11 O12 O13 ∑(O)

10.847 10.870 9.649 9.669 9.676 9.608 38.602

11.218 11.111 9.632 9.609 9.520 9.536 38.297

11.192 10.867 9.613 9.631 9.695 9.546 38.485

a

Data highlighted in bold refer to the most relevant changes upon deformation.

structure is deformed according to the two investigated modes. Though the observed variations in the bond type are substantial, there is likely no significant weakening of the structure because, when the bond length increases, the modification of the bond type from more polar to less polar tends to oppose decomposition. In short, the investigated deformations are certainly not responsible of any bond-rupture in the premelting T region. All calculations performed in this work suggest as premelting effects in sodium metasilicate (Na2SiO3) are not related to any “anomalous” vibrational mode in the high-frequency range of the substance. Very recently,56 a molecular dynamics simulation carried out on sodium disilicate (Na2Si2O5) showed that Na diffusivity strongly depends upon the sheet-like structure of the crystalline phase and leads in fact to premelting effects due to the enhanced dynamics of cations at T ≈ 800 K. The high diffusivity of Na atoms in the premelting range leads in turn to structural collapse of the Na subnetwork, whereas the silicate chains framework, although extensively deformed, is still preserved up to melting conditions. Owing to the structural similarities between Na2SiO3 and Na2Si2O5, it seems plausible that the same mechanism holds for both compounds, even though at slightly different temperature conditions. For Na2SiO3, the onset of premelting is observed 160 K below the room-pressure congruent melting point, when rearrangements of the Si2O6 chains begin to take place as seen in Raman and 29Si NMR spectra2,3,57 well after sodium mobility has itself set in at a much lower temperature of about 850 K. The bias between computed and calorimetric isobaric heat capacity at 1300 K is ∼5.9R, notably higher than the contribution arising eventually from translational freedom of Na atoms (∼3R) and difficult to explain when assuming the homogeneity of the crystal structure.

R = bond distance; ρbcp = electron density at bond critical point; δ(X,Y) = delocalization index. Data highlighted in bold refer to the changes upon deformation that are deemed to be more relevant for our study and in most cases also somewhat unexpected.

a

and vectorially deformed according to the ν54 and ν64 modes (Figure 8). The corresponding bond critical point (bcp) electron density and delocalization indexes values are also listed. In some Si−O bonds (bold character in Table 9), one observes significant variations in the electron density value induced by the deformation. For example, for the Si−O bond between centers 5 and 12 the bcp density value is reduced to half in ν54 and is tripled in ν64. The variation of electron density at the critical point and bond length are not systematically correlated as one would expect. When subjected to ν54 deformation, the Si−O bond length between centers 14 and 12 increases substantially but the delocalization index increases as well. The increase in bond length for the same bond when subject to the deformation mode ν64 is even higher, but the delocalization index is not apparently much affected. We may explain this intriguing behavior by assuming that the nature itself of the bond changes progressively as a result of the deformation: with the increase of the distance, the electron sharing (hence the degree of covalency) increases and less charge separation between Si and O occurs. Indeed, one may see (Table 10) that the atomic population of the involved atoms changes toward a lower ionicity (Si less positive, O less negative) when the

6. CONCLUSIONS Using the hybrid B3LYP density functional method, we computed the ab initio thermodynamic and thermophysical properties of Na2SiO3 in the Cmc21 structural state. The static properties were first evaluated through a symmetry-preserving relaxation procedure and then we computed, on the equilibrium structure, harmonic vibrational modes at the long-wavelength limit corresponding to the center of the Brillouin zone (k → 0). We obtained the acoustic modes from the elastic constant tensor by solving the Christoffel determinant. The acoustic branches were assumed to follow a sine-wave dispersion when traveling within the Brillouin zone. The procedure generated as well several thermophysical properties of interest in materials science and geophysics (transverse and longitudinal-wave velocities, shear modulus, Young modulus, Poisson ratio). All thermodynamic properties that depend on vibrational frequencies viz, heat capacities, thermal expansion, thermal derivative of the bulk M


Article


to D.B. (Visiting Professor Program 2016) is also warmly acknowledged.

modulus, thermal correction to internal energy, enthalpy, Gibbs free energy, were computed on the basis of quasi-harmonic mode-γ analysis of the volume effects on vibrational frequencies. Experiments and first-principles calculations are in substantial agreement regarding the Raman spectra of the Cmc21 phase at low T; however, Richet et al.2 noted that the Raman-active mode ν12 (two-components ν12 and ν12a) exhibits a nonlinear shift with temperatures above ∼900 K. They also observed that bands ν12 and ν17 (assigned respectively by the Authors to Si−O−Si bending + stretching and stretching motions at ν12 = 588 cm−1, ν17 = 965 cm−1) show a sudden increase in the area ratio with T about 100 K below the melting point. In the aim of investigating the isotropic chemical shifts of 23Na and 29Si in Na2SiO3 and to assess better the nature and implications of these two peculiar vibrations, we abstracted from the variational structure a cluster corresponding to an elongated packing along the c axis. The Ramsey contributions to the JNM coupling tensor between sodium and oxygen in Cmc21 cluster [Na8Si4O13]2− adopting the GIAO method are dominated by the interaction between the magnetic moments associated with the nuclear spins (i.e., the Fermi Contact FC). The JNM coupling tensor between sodium and oxygen is clearly affected by the interionic distance and decidedly more intense when nonbridging oxygen centers are involved. The calculated isotropic chemical shifts of 23Na is slightly higher than what was observed in the crystalline substance at room temperature3 (δiso ≈ 25) although substantially comparable. The calculated isotropic shift of 29Si is again in substantial agreement with what experimentally observed3 (δiso = −79.6). The spectral analysis carried out at the same level of theory adopted for the crystal assigns to the band ν54 of the cluster, which corresponds to the ν29 mode of the crystal, a frequency of 597 cm−1. The mode is effectively a mixed Si−O bending + Si− O−Si stretching motion affecting mainly the tetrameric chain. The ν64 mode of the cluster (954.2 cm−1) corresponds to the ν34 of the crystal (954.8 cm−1). This vibrational mode is essentially described by a stretching motion mostly affecting the bridging oxygens of the chain.2 To investigate further the role of deformation on structural stability, we submitted the cluster to an electron localization/delocalization analysis. In some Si−O bonds one observes significant variations in the electron density at the bond critical point induced by the deformation. The nature itself of the bond changes progressively as a result of the deformation: with the increase of the bond distance, the degree of bond covalency increases and the charge separation between Si and O lowers. Though the observed variations in the bond type are substantial, there is no significant weakening of the structure, because when the bond length increases, the modification of the bond type from more polar to less polar contrasts decomposition. In short, the investigated deformations are certainly not responsible of any bond-rupture in the premelting T region.

■

■

REFERENCES

(1) Richet, P.; Ingrin, J.; Mysen, B. O.; Courtial, P.; Gillet, P. Premelting effects in minerals: an experimental study. Earth Planet. Sci. Lett. 1994, 121, 589−600. (2) Richet, P.; Mysen, B. O.; Andrault, D. Melting and premelting of silicates: Raman spectroscopy and X-ray diffraction of Li2SiO3 and Na2SiO3. Phys. Chem. Miner. 1996, 23, 157−172. (3) George, A. M.; Richet, P.; Stebbins, J. F. Cation dynamics and premelting in lithium silicate (Li2SiO3) and sodium metasilicate (Na2SiO3): a high-temperature NMR study. Am. Mineral. 1998, 83, 1277−1284. (4) Grund, P. A.; Pizy, M. M. Structure cristalline du metasilicate de sodium anhydre, Na2SiO3. Acta Crystallogr. 1952, 5, 837−840. (5) McDonald, W. S.; Cruickshank, D. W. J. A reinvestigation of the structure of sodium metasilicate. Acta Crystallogr. 1967, 22, 37−43. (6) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648−5652. (7) Lee, C.; Yang, E.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (8) Peintinger, M. F.; Vilela Oliveira, D.; Bredow, T. Consistent Gaussian basis sets of triple-zeta valence with polarization quality for solid-state calculations. J. Comput. Chem. 2013, 34, 451−459. (9) Dovesi, R.; Roetti, C.; Freyria-Fava, C.; Prencipe, M. On the elastic properties of lithium, sodium and potassium oxide. An ab initio study. Chem. Phys. 1991, 156, 11−19. (10) Dovesi, R.; Saunders, V. R.; Roetti, C.; Orlando, R.; ZicovichWilson, C. M.; Pascale, F.; Civalleri, B.; Doll, K.; Harrison, N. M.; Bush, I. J.; et al. CRYSTAL14 User’s Manual; University of Torino: Torino, Italy, 2014. (11) Kohanoff, J. Electronic Structure Calculations for Solids and Molecules: Theory and Computational Methods; Cambridge University Press: Cambridge, U.K., 2006. (12) Ottonello, G.; Civalleri, B.; Ganguly, J.; Vetuschi Zuccolini, M.; Noël, Y. Thermophysical properties of the α-β-γ polymorphs of Mg2SiO4: a computational study. Phys. Chem. Miner. 2009, 36, 87−106. (13) Lin, C. C.; Leung, K. S.; Shen, P.; Chen, S. F. Elasticity and structure of the compounds in the wollastonite (CaSiO3)-Na2SiO3 system: from amorphous to crystalline state. J. Mater. Sci.: Mater. Med. 2015, 26, 39. (14) Doll, K. Implementation of analytical Hartree-Fock gradients for periodic systems. Comput. Phys. Commun. 2001, 137, 74−88. (15) Pascale, F.; Zicovich-Wilson, C. M.; Lopez-Gejo, F.; Civalleri, B.; Orlando, R.; Dovesi, R. The calculation of the vibrational frequencies of crystalline compounds and its implementation in the CRYSTAL code. J. Comput. Chem. 2004, 25, 888−897. (16) De La Pierre, M.; Belmonte, D. Ab initio investigation of majorite and pyrope garnets: lattice dynamics and vibrational spectra. Am. Mineral. 2016, 101, 162−174. (17) Barrow, G. M. Introduction to Molecular Spectroscopy; McGrawHill: New York, 1962. (18) Hess, B. A.; Schaad, L. J.; Carsky, P.; Zahradnik, R. Ab initio calculations of vibrational spectra and their use in the identification of unusual molecules. Chem. Rev. 1986, 86, 709−730. (19) Dall’Olio, S.; Dovesi, R.; Resta, R. Spontaneous polarization as a Berry phase of the Hartree-Fock wave function: The case of KNbO3. Phys. Rev. B: Condens. Matter Mater. Phys. 1997, 56, 10105−10114. (20) Noël, Y.; Zicovich-Wilson, C. M.; Civalleri, B.; D’Arco, Ph.; Dovesi, R. Polarization properties of ZnO and BeO: an ab initio study through the Berry phase and Wannier functions approaches. Phys. Rev. B: Condens. Matter Mater. Phys. 2001, 65, 014111. (21) Maschio, L.; Kirtman, B.; Rérat, M.; Orlando, R.; Dovesi, R. Ab initio analytical Raman intensities for periodic systems through a coupled perturbed Hartree-Fock/Kohn-Sham method in an atomic orbital basis. I. Theory. J. Chem. Phys. 2013, 139, 164101.

AUTHOR INFORMATION

Corresponding Author

*Giulio Ottonello. E-mail: [email protected]. Phone: (0039) 353 010 8136. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS G.O. is indebted with IPGP (Paris) for kind hospitality. Financial support by the Institut de Physique du Globe de Paris (IPGP) N


Article


N.; et al. Gaussian 03, Revision B.05; Gaussian, Inc.: Wallingford, CT, 2004. (44) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab initio calculation of vibrational absorption and circular dichroism spectra using density functional force fields. J. Phys. Chem. 1994, 98, 11623−11627. (45) Tomasi, J.; Persico, M. Molecular interactions in solution: An overview of methods based on continuous distributions of the solvent. Chem. Rev. 1994, 94, 2027−2094. (46) Gatti, C.; Ottonello, G.; Richet, P. Energetics and bonding in aluminosilicate rings with alkali metal and alkaline-earth metal chargecompensating cations. J. Phys. Chem. A 2012, 116, 8584−8598. (47) Vayssilov, G. N.; Rösch, N. Density functional studies of alkaliexchanged zeolites: Basicity and core-level shifts of framework oxygen atoms. J. Catal. 1999, 186, 423−432. (48) Vayssilov, G. N.; Lercher, J. A.; Rösch, N. Interaction of methanol with alkali metal exchanged molecular sieves. 2. Density functional study. J. Phys. Chem. B 2000, 104, 8614−8623. (49) Vayssilov, G. N.; Staufer, M.; Belling, T.; Neyman, K. M.; Knozinger, H.; Rösch, N. Density functional studies of alkali-exchanged zeolites. Cation location at six-rings of different aluminum content. J. Phys. Chem. B 1999, 103, 7920−7928. (50) Bader, R. F. W.; Stephens, M. E. Fluctuation and correlation of electrons in molecular systems. Chem. Phys. Lett. 1974, 26, 445−449. (51) Poater, J.; Solà, M.; Duran, M.; Fradera, X. The calculation of electron localization and delocalization indices at the Hartree-Fock, density functional and post-Hartree-Fock levels of theory. Theor. Chem. Acc. 2002, 107, 362−371. (52) Matito, E.; Solà, M.; Salvador, P.; Duran, M. Electron sharing indexes at the correlated level. Application to aromaticity calculations. Faraday Discuss. 2007, 135, 325−345. (53) Fradera, X.; Austen, M. A.; Bader, R. F. W. The Lewis model and beyond. J. Phys. Chem. A 1999, 103, 304−314. (54) Poater, J.; Duran, M.; Solà, M.; Silvi, B. Theoretical evaluation of electron delocalization in aromatic molecules by means of Atoms in Molecules (AIM) and Electron Localization Function (ELF) topological approaches. Chem. Rev. 2005, 105, 3911−3947. (55) Gatti, C.; Lasi, D. Source function description of metal-metal bonding in d-bond organometallic compounds. Faraday Discuss. 2007, 135, 55−78. (56) Mantisi, B.; Micoulaut, M. Premelting and cation mobility in simple silicates: contrasting dynamics in the crystalline and molten state. J. Non-Cryst. Solids 2016, 440, 1−6. (57) Nesbitt, H. W.; Bancroft, G. M.; Henderson, G. S.; Richet, P.; O’Shaughnessy, C. Melting, crystallization and the glass transition: Toward a unified description for silicate phase transitions. Am. Mineral. 2016, in press, DOI: 10.2138/am-2017-5852. (58) Dovesi, R.; Ermondi, E.; Ferrero, E.; Pisani, C.; Roetti, C. HartreeFock study of lithium hydride with the use of a polarizable basis set. Phys. Rev. B: Condens. Matter Mater. Phys. 1984, 29, 3591−3600.

(22) Maschio, L.; Kirtman, B.; Rérat, M.; Orlando, R.; Dovesi, R. Ab initio analytical Raman intensities for periodic systems through a coupled perturbed Hartree-Fock/Kohn-Sham method in an atomic orbital basis. II. Validation and comparison with experiments. J. Chem. Phys. 2013, 139, 164102. (23) Born, M.; Huang, K. Dynamical Theory of Crystal Lattices; Oxford University Press: Oxford, U.K., 1954. (24) Musgrave, M. J. P. Crystal Acoustics; Holdan-Day: Boca Raton, FL, 1970. (25) Nye, J. F. Physical Properties of Crystals; Oxford University Press: Oxford, U.K., 1957. (26) Mainprice, D. In Treatise on Geophysics; Price, G. D., Ed.; Elsevier: Amsterdam, 2007; Vol. 2, pp 437−491. (27) Kieffer, S. W. Thermodynamics and lattice vibrations of minerals: 1. Mineral heat capacities and their relationships to simple lattice vibrational models. Rev. Geophys. 1979, 17, 1−19. (28) Hill, R. W. The elastic behaviour of a crystalline aggregate. Proc. Phys. Soc., London, Sect. A 1952, 65, 349−354. (29) Erba, A.; Mahmoud, A.; Belmonte, D.; Dovesi, R. High pressure elastic properties of minerals from ab initio simulations: the case of pyrope, grossular and andradite silicate garnets. J. Chem. Phys. 2014, 140, 124703. (30) Mainprice, D. A Fortran program to calculate seismic anisotropy from the lattice preferred orientation of minerals. Comput. Geosci. 1990, 16, 385−393. (31) Ottonello, G.; Vetuschi Zuccolini, M.; Civalleri, B. Thermochemical and thermo-physical properties of stishovite: An ab-initio allelectron investigation. CALPHAD: Comput. Coupling Phase Diagrams Thermochem. 2009, 33, 457−468. (32) Ottonello, G.; Civalleri, B.; Ganguly, J.; Perger, W. F.; Belmonte, D.; Vetuschi Zuccolini, M. Thermo-chemical and thermo-physical properties of the high-pressure phase anhydrous B (Mg14Si5O24): An ab-initio all-electron investigation. Am. Mineral. 2010, 95, 563−573. (33) Ottonello, G.; Vetuschi Zuccolini, M.; Belmonte, D. The vibrational behaviour of silica clusters at the glass transition: Ab initio calculations and thermodynamic implications. J. Chem. Phys. 2010, 133, 104508. (34) Belmonte, D.; Ottonello, G.; Vetuschi Zuccolini, M. Melting of α-Al2O3 and vitrification of the undercooled alumina liquid: Ab initio vibrational calculations and their thermodynamic implications. J. Chem. Phys. 2013, 138, 064507. (35) Belmonte, D.; Ottonello, G.; Vetuschi Zuccolini, M. Ab initio thermodynamic and thermophysical properties of sapphirine endmembers in the join Mg4Al8Si2O20-Mg3Al10SiO20. Am. Mineral. 2014, 99, 1449−1461. (36) Kieffer, S. W. Thermodynamics and lattice vibrations of minerals: 2. Vibrational characteristic of silicates. Rev. Geophys. 1979, 17, 20−34. (37) Kieffer, S. W. Thermodynamics and lattice vibrations of minerals: 3. Lattice dynamics and an approximation for minerals with application to simple substances and framework silicates. Rev. Geophys. 1979, 17, 35−59. (38) Ottonello, G.; Attene, M.; Ameglio, D.; Belmonte, D.; Vetuschi Zuccolini, M.; Natali, M. Thermodynamic investigation of the CaOAl2O3-SiO2 system at high P and T through polymer chemistry and convex-hull techniques. Chem. Geol. 2013, 346, 81−92. (39) Kelley, K. K. The specific heats at low temperatures of crystalline ortho-, meta-, and di-silicates of sodium. J. Am. Chem. Soc. 1939, 61, 471−473. (40) Barin, I. Thermochemical Data of Pure Substances; VCH: Weinheim, 1995. (41) Chase, M. W., Jr. NIST − JANAF Thermochemical Tables. Journal of Physical and Chemical Reference Data, Monograph No. 9; American Chemical Society and American Institute of Physics for the National Institute of Standards and Technology: New York, 1998. (42) Richet, P.; Bottinga, Y.; Téqui, C. Heat capacity of sodium silicate liquids. J. Am. Ceram. Soc. 1984, 67, C6−C8. (43) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. O


Ab Initio Thermodynamic and Thermophysical Properties of Sodium

Recommend Documents