- Email: [email protected]
Quantitative interpretation of NMR spectra of disordered solids
SOLID STATE
Solid State Ionics 53-56 (1992) 868-~74 North-Holland
IOHICS Quantitative interpretation of NMR spectra of disordered solids Pacifico Cofrancesco
1, Sergio Scotti, Marco Villa
Dipartimento di Fisica "A. Volta". GNSM-CISM, Via Bassi 6, 27100 Pavia, Italy
Piercarlo Mustarelli Dipartimento di Chimica-Fisica, Viale Taramelli 16, 27100 Pavia, ltalv
and Ales Gottvald Institute of Scientific Instruments, Czechoslovak Academy of Sciences, Krcilovopolskd 147, 61264 Brno, (2echoslovakia
Profile analysis of quadrupole-broadened NMR spectra of disordered solids is a complex and computationally intensive task that we have undertaken by developing general purpose optimization programs. After presenting the spectroscopic problem (firstor second-order quadrupolar effects) and the different types of optimizers (deterministic or stochastic), we discuss how efficiently the "'solution" of a simple problem can be found, how a model of disorder can be introduced and progressively refined, which kind of information can be obtained. As example, ~B magic angle spinning NMR spectra in vitreous borates are quantitatively interpreted.
1. Introduction Over two-thirds o f the stable nuclear isotopes have spin I > ½, and carry an electric quadrupole m o m e n t Q. The coupling o f Q with electric field gradients ( E F G s ) at the nuclear position affects the dynamics of the magnetic m o m e n t associated with the angular m o m e n t u m hi, and can be studied with nuclear magnetic resonance ( N M R ) techniques to yield inform a t i o n about coordination and dynamics. The local nature of this interaction makes it a powerful probe in structural investigations o f microcrystalline and disordered materials, where usefulness o f diffraction techniques is limited. However, relative little use has so far been m a d e of the q u a d r u p o l a r probe, the m a i n reasons being the following: ( a ) The q u a d r u p o l a r interaction averages out in liquids, while in solids it m a y be obscured by other interactions, or it may yield spectra too wide to be observed in full. Author to whom correspondence should be addressed.
( b ) When the EFGs are known, calculation o f a powder pattern is tedious but straightforward; solution o f the inverse problem is either an impossible or a very complex task. The N M R spectroscopists prefer studies of spin I = ½ nuclei [1,2], where the chemical shift interaction can easily be studied: spectral width is seldom a problem, and we usually analyze only the rotationally invariant (isotropic) part of the chemical shift, which enormously simplifies spectral assignments or calculations. We argue that this situation is going to change, as better experimental techniques overcome sensitivity problems, and computational tools are developed, which ease a quantitative interpretation of complex spectra. In particular, all q u a d r u p o l a r nuclei but seven have a semi-integer spin: their central transition ( ½,--,- ½) in a strong external field B shows second-order q u a d r u p o l a r broadening which scale with 1 / B , and can be further reduced by spinning the sample at the "magic angle" ( M A S ) . This rotation averages out the dipolar and the anisotropic chem-
0167-2738/92/$ 05.00 © 1992 Elsevier Science Publishers B.V. All rights reserved.
P. Cofrancesco et al. / Quantitative interpretation of NMR spectra of disordered solids
ical shift interactions, and, under suitable experimental conditions, yields a spectrum dominated by quadrupolar effects. A privilege of the N M R spectroscopy is that we can build mathematical models of the spectrum directly from the first principles, while introducing essentially no approximations. Then, we have to optimize the choice of the adjustable parameters of a model, and compare different models. A range of optimization techniques are currently being suggested or developed which will have a profound influence upon data interpretation and engineering. An overview of the field will be given in section 2, where we introduce the categories of the deterministic and stochastic optimizers. In section 3 we study how efficiently these two different types of optimizers succeed in retrieving the EFG parameters from a simulated MAS spectrum broadened by second-order quadrupolar effects. Then, two actual ~~B spectra in borate glasses are discussed to show an example of model-refining and to comment upon the information provided by our analysis.
2. Deterministic and stochastic optimizers Let F~(co) be an experimental function and Fc(o), xb x2 ..... x,) a calculated function, which depends upon the vector of parameters x = (xi) (degrees of freedom ). The choices of the xi may be limited (e.g., by conditions such as a,
(i )
Usually, F~ is known over a set of points {co,), and the norm in eq. ( 1 ) is computed as [F~ (09,) -F~(~o,, x ) 1 2 .
(2)
i
Many optimizers build a sequence o f " s o l u t i o n s " x (k) with the recursive equation x(k+ J) = X ( k ) + a (k)s (k) ,
(3)
beginning with an initial guess x (°). Here k is the number of iteration, a Ck) is called the step, and s (k) [he "slope". Different choices of the slope vector
869
characterize the various methods; for example, if we move our point in the direction of the fastest decrease o f f s(k~= - g r a d f ( x (k)) ,
(4)
we have the so-called steepest-descent method. All popular choices of the slope vector (e.g. quasi-Newton, conjugate gradient, etc) involve a finite-difference estimation of the derivatives of the error function in x (k). All these methods may be grouped together in the large class of the so-called higher-order deterministic optimization methods (HOD O M ) . We have developed [3] a H O D O M which, at each iteration, attempts to find a new point x' by moving in the direction of the negative gradient by a step not related to the gradient size, using 1D linesearch. The new point should belong to the allowed domain of x, and the displacement I x ' - x ~k) ] should be smaller than an iteration-dependent limit l (k). If f ( x' ) < f ( x ( k ) ) , x' becomes the new trajectory point x (~+ ~), and a new iteration begins with [(~+~) slightly larger than l (h~. Otherwise, the limit l (k~ is decreased (by a small fraction), and the search around x (k) continue until the maximum step size, a n d / o r the difference I f ( x ' ) - f ( x ~k~ )l, fall below prescribed values. In a constrained optimization problem, the x (k) trajectory initially explores the boundaries of the parameter's domain, and can easily escape from local minima o f f ( x ) if the percentage adjustment of the step limit is made very small. Obviously, the smaller the change of the move limit, the slower the convergence to a solution. Furthermore, computation of the gradient needs n + 1 function calls, where n is the number of degrees of freedom. In the optimization of the signal shape, a new class of stochastic methods may be usefully introduced. These methods are based on a technique known as s i m u l a t e d annealing [4,5]; and even if they have c o m m o n mathematical foundations [8], they have different names (genetic algorithms [6], evolution strategies [7], ...) and different physical interpretations for the different research areas in which they are used. We generally can call them zeroth-order stochastic optimization methods (ZOSOM). In such methods, we start the kth iteration with a set of/z "parents" vectors: x~ ~ ..... x(,k~, and generate 2 random "children" vectors y~"), ..., y~k) constrained to a probability density function p. For example, the
870
P. Cofrancesco et al. /Quantitative interpretation oI'NMR spectra o/'disordered solids
probability to have a child near point x = (xi) may be chosen p r o p o r t i o n a l to the m u l t i d i m e n s i o n a l Gaussian probability density function
(
p(x,m,d)=exp
- ~
1,
(5)
where m, is a mean over the parents x , and d~ are the components o f a dispersion vector adjusted at each iteration. We have a (/2+2) strategy if in the next iteration the new parents are the # vectors with the smallest errors in the set comprising both parents and children. We have a (#, 2) strategy if the best # vectors are selected among the " c h i l d r e n " only ( 2 > #). The so-called Monte Carlo m e t h o d is an example of a ( 1 + 1 ) strategy. Fig. 1 graphically illustrates the (3 + 4) and the (2, 6) strategies. The other ingredient o f an evolution strategy is the rule by which the dispersions are modified; if the parents are close to the m i n i m u m and the range of search is much larger than the distance from it, only a small fraction of children will be better than the parents, and Idl should be reduced; when children are frequently better than parents the dispersion should be increased, to explore regions farther away from m. It follows that a decision about the dispersion at the next step should be taken on the basis o f a statistically significant n u m b e r of generations. 3
~=4
,L
£
g
p
X
X
X
- - >
(2)
Y
_
__
Y
Y
Y
(a) ~t 2
L=6 p
(I)
x
x
> Y
Y (.J) Y
YjLY)
(2)
X
X
> Y
Y
Y
Y
Y
Y
(b) Fig. 1. Graphic illustration of the ( 3 + 4) and the (2, 6 ) strategies.
Tuning elements of an evolution strategy are, typically, the length of the evolution histoD,, the fractional change of dispersion, the critical success rate o f children which triggers an increase of the search range. With increasing the dimension o f the x-domain we are likely to need to increase the evolution history,, but total c o m p u t a t i o n time may increase much less. The r a n d o m search explores simultaneously all x-dimensions, and the o p t i m i z e r overhead is almost independent of the n u m b e r of parameters. These stochastic methods have many attractive features, when c o m p a r e d with deterministic ones. Firstly, they are more effective in reaching a global minimizer, or at least a very stable local minimizer: secondly, the accuracy of solution is high because it does not d e p e n d upon the accuracy of the gradients: thirdly, the method can be applied even if the derivative of the objective function does not exist: and finally the treatment of constraint conditions is usually much simpler than in HODOMs. From the timeefficiency point of view, some results show that stochastic methods are comparable with deterministic ones, and sometimes they work even better [8,9].
3. T e s t i n g t h e o p t i m i z e r s
No optimizer warrants the convergence to the absolute m i n i m u m o f f ( x ) , nor every' problem can be treated with an optimizer. It has been argued [ 10] that d e t e r m i n a t i o n o f the E F G from a second-order quadrupolar spectrum is one of these problems, since substantially different EFGs may produce very similar powder patterns. An " e l e m e n t a r y " N M R MAS profile can be calculated at a given field B by assigning the strength of the q u a d r u p o l a r interaction VQ, the asymmetry p a r a m e t e r t/, and the width A of a broadening function (a Gaussian in our case). We have tested the optimizers by using as "experimental" functions spectra c o m p u t e d for known values of VQ, ~1, A. By trial and error we found the move limit adj u s t m e n t factor (0.95) which ensures convergence of our H O D O M to the "right" parameters, no matter how the initial point is chosen. In the process we identified several of the potential ambiguities denounced by Snyder et al. [ 10 ]. However, there is no doubt that with noiseless and undistorted data a ro-
P. Cofrancesco et aL /Quantitative interpretation of NMR spectra of disordered solids
871
Table 1 Data on optimization runs with a 1 + 1 evolution strategy q
F-calls
Errors of solution
18ql
8~,Q(%)
0.00 0.00
0.120 0.034
0.154 0.000
0.03 0.16
0.24 0.05
109 270
0.33 0.33
0.044 0.053
0.380 0.310
10.80 14.00
0.67 0.88
135 127
0.66 O.66
0.010 O.O21
0.460 1.230
5.04 1.36
0.46 1.15
400 147
1.00 1.00
0.082 0.053
2.380 1.530
3.72 1.67
0.35 0.22
167 269
bust optimizer will eventually approach the true solution. A typical solution has a rms< 0.1% and requires about 103 function calls. Data about optimization runs with a 1 + 1 evolution strategy are summarized in table 1. Each optimization has been repeated two times for spectra computed at q=0, 0.33, 0.66 and 1. The table evidences the order-of-magnitude fluctuations of the accuracy of the solution expressed as rms deviation, and over a 300% variation in the number of function calls. This fact had to be expected, since both starting point and evolution depend upon the rolling of a dice, and the random number generator is casually seeded before each run. Runs with about 100 function calls are typical, but accuracy in determination of the parameters (in particular q) is not always satisfactory. One may tune the parameters of the optimizer to reach better solutions in longer runs; alternatively, one may repeat several times the same run to estimate experimentally how accurate the solutions are, and what is the sensitivity of the error function near the minimum to the various parameters.
4. Two examples With a robust and reliable optimizer at our disposal, our attention should focus upon building an appropriate model of the experimental spectrum. In a polycrystalline powder with the resonant nuclei all in equivalent positions we may use the single-site (UQ, r/, A)-model used above. A should be strongly con-
8A (%)
rms (%)
strained with independent estimates of the instrumental broadening, and all causes of distortions discounted from the experimental data, or inserted in the model. To simulate the glassy counterpart of this solid we have to compromise between a realistic picture of disorder and the simplicity of the mathematical description. We may begin with the simplest of the models of disorder, and progressively refine it until a "reasonable" agreement with the experiment has been reached. An example of this procedure is illustrated for the case of the ~IB MAS spectrum of vitreous B203, chosen because it was expected to consist only of borons bonded to three bridging oxygens (B3). In model A we have assumed that disorder just results in an additional broadening, of the same nature of the instrumental one, and lifted the constraint upon A. Our HODOM gives a A about three times the instrumental one, but the rms error is high, 9.09%. Model B assumes that disorder is described by a normal distribution of q values, constrained between 0 and 1, centered at qo and with a width 2At/ at l / e P(qo). However, a substantially better agreement is reached (model C, fig. 2) by normally distributing the quadrupolar coupling constant PQ in a way statistically independent of the q-distribution:
P("Q)
L \ ~z,Q ] =]A
(6)
The broadening reaches the more reasonable value of 164 Hz, and the error reduces to 5.06%. However,
872
P. Cofrancesco et al. /Quantitative interpretation oJ'NMR spectra of disordered solids
tO moisture contamination through the reaction 2 B 3 + H20,--.2 ( B 4 - + H + ) ,
2,83
KHz
-2.88
Fig. 2. Model C. the difference plot Fe(Og) - F c ( o g ) (fig. 2) clearly indicates that the experimental spectrum contains small and narrow components, which we model with a Lorentzian-Gaussian profile 1
{l-g[
( ~ ) ]
2
--1
2
+0.69 g e x p [ - 0 . 6 9 ( ~o-)--- )e)o
]},
(7)
and four additional parameters: p, relative amplitude of the narrow component; ~Oo, position of the baricentre; g, relative weight of the Gaussian component; A', half height at half width of L + G functions. The parameters relative to the B3 sites remain essentially equal to those of model C, but the rms error reduces to 2.41% (model D). Position and width of the narrow component agree with the values characteristic of the component arising from fourcoordinated borons (B4). If this component is due
its amplitude ( p = 3%) allows to evaluate a reasonable water content of 0.7 wt%. Table 2 lists optimum values of parameters for models A, B, C and D. Fig. 3 shows an ~B MAS spectrum (open circles) obtained in the Li + glassy conductor Li20:B203: Li2SQ. The modelled one (solid curve) has been obtained with our H O D O M by assuming that four different types of structural units are present. (i) Symrnetric B3 units, with the same profile of the B3 units in B203. The only adjustable parameters of this profile are its weight p(B3S) and chemical shift. (ii) B4 6nits, with a small quadrupolar broadening, and a profile described by eq. (7). Beside the weight factor p(B4), the adjustable parameters are position (~o0) , width (A'), and fraction (g) of the Gaussian character of the profile. (iii) Asymmetric B3 units, with borons connected to one or two non-bridging oxygens. In an extensive N M R study of lithium-rich borate glasses, Yun and Bray [ l 1 ] have been able to determine average values Of UQand ~/for this unit; we have derived its EFG parameters by fitting our MAS spectra o f L i 2 0 : B 2 0 3 glasses. We obtained UQo= 1.38 MHz, and qo=0.43, values consistent with the values quoted by Yun and Bray [ 11 ]. In addition, we found a very small dispersion for vo ( A u o = 3 0 kHz), and a very large one for r/(At/= 1.96). As for the (i) component, the only adjustable parameters for the B3A profile are weight, p ( B 3 A ) , and chemical shift. (iv) Borate units bonded to sulphur, SB, with weight p(SB) = 1 - p ( B 3 S ) - p ( B 4 ) - p ( B 3 A ) . Their possible structures have been discussed in a previous
Table 2 Optimum values of the parameters for A, B, C and D models, Model
qo
Aq
UQo(MHz)
AUQ
A (kHz)
rms (%)
A B C
0.22 0 0
1.799 1.860
1.357 1.357 1.357
0.043
0.327 0.278 0.164
9.09 7.99 5.06
D
coo (kHz)
g
A' (kHz)
p(B4)
rms (%)
0.808
0.89
0.337
0.03
2.41
P. Cofrancesco et al. /Quantitative interpretation of NMR spectra of disordered solids
6.71
I(]-Iz
2.16
Fig. 3. ~B MAS spectrum (see text). Table 3 Optimal parameters p (B3S)=0.151 p (B3A)=0.145 p (B4) =0.416 090=4.908 kHz 4' =0.208 kHz g=0.659
p (SB) =0.288 t/o= 0.6 A~/= 2.89 UQo= 1.148 MHz AvQ= 1.52 zi=0.194 kHz
work [12]. Here, we estimate the EFG parameters of SB by letting the optimizer search for the best fit values of VOo, Az,o, qo and Aq, with essentially no constraints. With this model, a good agreement with experiment, and a low rms error of 3.54% is reached. Table 3 lists all data obtained from this optimization run, with the exclusion of the chemical shift parameters. We acknowledge that this model requires the checking of several assumptions (not discussed here) and implies the knowledge of complex profiles obtained from independent measurements; however, all this effort is justified by the rich set of structure-related parameters we have been able to extract from a single spectrum.
6. Conclusions A primary goal of the paper is that of conveying the idea that profile analysis of quadrupolar spectra is a daunting but rewarding task. With an ideal ex-
873
perimental setup we may plan MAS experiments at very high polarizing fields, which allow estimating number and amplitudes of the distinct components by making the chemical shift interaction predominant, and at relatively low magnetic fields, where we are more sensitive to the EGF distribution of the individual components. However, the second-order quadrupolar interaction is always associated with an overall displacement of the line, which we handle by simultaneously optimizing EFG and chemical shift parameters. Consistency of the latter ones over families of related glasses is a strong evidence in favour of the model chosen. When it can be proven that a family of glasses consists of a small number of structural units with characteristic EFG distributions, the analysis simplifies enormously, since we can adopt as spectral components the characteristic profiles of these units, and adjust their amplitudes to fit the experimental data. In general, however, a model results from linear combination of many primitive profiles, each representing the powder pattern corresponding to a choice of r/, u, and a. The computation time may become an issue even with the most powerful hardware currently available. We deal with it by selecting and tuning the optimization technique best suited to a given problem, and by writing a most efficient code for the error function. In N M R spectral modeling, the technique used to compute this function plays a very critical role. For this reason, we have evaluated several general methods for calculating, with prescribed resolution and accuracy, the "powder averages" leading to a primitive profile. The results of these studies will be presented elsewhere.
References [ 1 ] G. Engelhardt and D. Michel, High Resolution Solid State NMR of Silicates and Zeolites (Wiley, New York, 1987). [2 ] G.L. Turner, R.J. Kirkpatrick, S.H. Risbud and E. Oldfield, Ceram. Bull. 66 (1987) 656. [3] P. Mustarelli, M. Rudnicki, A. Savini, F. Savoldi and M. Villa, Mag. Res. Imag. 8 (1990) 101. [4] S. Kirkpatrick, C.D. Gelatt and M.P. Vecchi, Science 220 (1983) 671. [5] I.O. Bohachevsky, M.E. Johnson and M.L. Stein, Technometrics 28 (1986) 209. [6] L. Davis, Genetic Algorithms and Simulated Annealing (Morgan-Kauffmann, Los Altos, 1989 ).
874
P. Cofrancesco et al. / Quantitative interpretation qf NMR spectra qf disordered solids
[7] A. Gottvald, K. Preis and M. Friedrich, OptimiZation Methods for Computational Electromagnetics, IGTE Syrup., Graz (1990) p. 23. [8] A. Gottvald, K. Preis, C. Magele, O. Biro and A. Savini, Global Optimization Methods for Computational Electromagnetics, in: Proc. Compumag, Sorrento, Italy, Vol. 2 (1991) pp. 1077-1080.
[9] J.H. Kalivas, N. Roberts and J.M. Sutter, Anal. Chem. 61 (1989) 2024. [10] L.C. Snyder, G.E. Peterson and C.R. Kurkjian. J. Chem. Phys. 64 (1976) 1569. [ 11 ] Y.H. Yun and P.J. Bray, J. Non-Cryst. Solids 44 ( 1981 ) 227. [ 12 ] P. Mustarelli, S. Scotti, M. Villa and P.R. Gandhi, Solid Stale lonics 39 (1990) 217.
Solid State Ionics 53-56 (1992) 868-~74 North-Holland
IOHICS Quantitative interpretation of NMR spectra of disordered solids Pacifico Cofrancesco
1, Sergio Scotti, Marco Villa
Dipartimento di Fisica "A. Volta". GNSM-CISM, Via Bassi 6, 27100 Pavia, Italy
Piercarlo Mustarelli Dipartimento di Chimica-Fisica, Viale Taramelli 16, 27100 Pavia, ltalv
and Ales Gottvald Institute of Scientific Instruments, Czechoslovak Academy of Sciences, Krcilovopolskd 147, 61264 Brno, (2echoslovakia
Profile analysis of quadrupole-broadened NMR spectra of disordered solids is a complex and computationally intensive task that we have undertaken by developing general purpose optimization programs. After presenting the spectroscopic problem (firstor second-order quadrupolar effects) and the different types of optimizers (deterministic or stochastic), we discuss how efficiently the "'solution" of a simple problem can be found, how a model of disorder can be introduced and progressively refined, which kind of information can be obtained. As example, ~B magic angle spinning NMR spectra in vitreous borates are quantitatively interpreted.
1. Introduction Over two-thirds o f the stable nuclear isotopes have spin I > ½, and carry an electric quadrupole m o m e n t Q. The coupling o f Q with electric field gradients ( E F G s ) at the nuclear position affects the dynamics of the magnetic m o m e n t associated with the angular m o m e n t u m hi, and can be studied with nuclear magnetic resonance ( N M R ) techniques to yield inform a t i o n about coordination and dynamics. The local nature of this interaction makes it a powerful probe in structural investigations o f microcrystalline and disordered materials, where usefulness o f diffraction techniques is limited. However, relative little use has so far been m a d e of the q u a d r u p o l a r probe, the m a i n reasons being the following: ( a ) The q u a d r u p o l a r interaction averages out in liquids, while in solids it m a y be obscured by other interactions, or it may yield spectra too wide to be observed in full. Author to whom correspondence should be addressed.
( b ) When the EFGs are known, calculation o f a powder pattern is tedious but straightforward; solution o f the inverse problem is either an impossible or a very complex task. The N M R spectroscopists prefer studies of spin I = ½ nuclei [1,2], where the chemical shift interaction can easily be studied: spectral width is seldom a problem, and we usually analyze only the rotationally invariant (isotropic) part of the chemical shift, which enormously simplifies spectral assignments or calculations. We argue that this situation is going to change, as better experimental techniques overcome sensitivity problems, and computational tools are developed, which ease a quantitative interpretation of complex spectra. In particular, all q u a d r u p o l a r nuclei but seven have a semi-integer spin: their central transition ( ½,--,- ½) in a strong external field B shows second-order q u a d r u p o l a r broadening which scale with 1 / B , and can be further reduced by spinning the sample at the "magic angle" ( M A S ) . This rotation averages out the dipolar and the anisotropic chem-
0167-2738/92/$ 05.00 © 1992 Elsevier Science Publishers B.V. All rights reserved.
P. Cofrancesco et al. / Quantitative interpretation of NMR spectra of disordered solids
ical shift interactions, and, under suitable experimental conditions, yields a spectrum dominated by quadrupolar effects. A privilege of the N M R spectroscopy is that we can build mathematical models of the spectrum directly from the first principles, while introducing essentially no approximations. Then, we have to optimize the choice of the adjustable parameters of a model, and compare different models. A range of optimization techniques are currently being suggested or developed which will have a profound influence upon data interpretation and engineering. An overview of the field will be given in section 2, where we introduce the categories of the deterministic and stochastic optimizers. In section 3 we study how efficiently these two different types of optimizers succeed in retrieving the EFG parameters from a simulated MAS spectrum broadened by second-order quadrupolar effects. Then, two actual ~~B spectra in borate glasses are discussed to show an example of model-refining and to comment upon the information provided by our analysis.
2. Deterministic and stochastic optimizers Let F~(co) be an experimental function and Fc(o), xb x2 ..... x,) a calculated function, which depends upon the vector of parameters x = (xi) (degrees of freedom ). The choices of the xi may be limited (e.g., by conditions such as a,
(i )
Usually, F~ is known over a set of points {co,), and the norm in eq. ( 1 ) is computed as [F~ (09,) -F~(~o,, x ) 1 2 .
(2)
i
Many optimizers build a sequence o f " s o l u t i o n s " x (k) with the recursive equation x(k+ J) = X ( k ) + a (k)s (k) ,
(3)
beginning with an initial guess x (°). Here k is the number of iteration, a Ck) is called the step, and s (k) [he "slope". Different choices of the slope vector
869
characterize the various methods; for example, if we move our point in the direction of the fastest decrease o f f s(k~= - g r a d f ( x (k)) ,
(4)
we have the so-called steepest-descent method. All popular choices of the slope vector (e.g. quasi-Newton, conjugate gradient, etc) involve a finite-difference estimation of the derivatives of the error function in x (k). All these methods may be grouped together in the large class of the so-called higher-order deterministic optimization methods (HOD O M ) . We have developed [3] a H O D O M which, at each iteration, attempts to find a new point x' by moving in the direction of the negative gradient by a step not related to the gradient size, using 1D linesearch. The new point should belong to the allowed domain of x, and the displacement I x ' - x ~k) ] should be smaller than an iteration-dependent limit l (k). If f ( x' ) < f ( x ( k ) ) , x' becomes the new trajectory point x (~+ ~), and a new iteration begins with [(~+~) slightly larger than l (h~. Otherwise, the limit l (k~ is decreased (by a small fraction), and the search around x (k) continue until the maximum step size, a n d / o r the difference I f ( x ' ) - f ( x ~k~ )l, fall below prescribed values. In a constrained optimization problem, the x (k) trajectory initially explores the boundaries of the parameter's domain, and can easily escape from local minima o f f ( x ) if the percentage adjustment of the step limit is made very small. Obviously, the smaller the change of the move limit, the slower the convergence to a solution. Furthermore, computation of the gradient needs n + 1 function calls, where n is the number of degrees of freedom. In the optimization of the signal shape, a new class of stochastic methods may be usefully introduced. These methods are based on a technique known as s i m u l a t e d annealing [4,5]; and even if they have c o m m o n mathematical foundations [8], they have different names (genetic algorithms [6], evolution strategies [7], ...) and different physical interpretations for the different research areas in which they are used. We generally can call them zeroth-order stochastic optimization methods (ZOSOM). In such methods, we start the kth iteration with a set of/z "parents" vectors: x~ ~ ..... x(,k~, and generate 2 random "children" vectors y~"), ..., y~k) constrained to a probability density function p. For example, the
870
P. Cofrancesco et al. /Quantitative interpretation oI'NMR spectra o/'disordered solids
probability to have a child near point x = (xi) may be chosen p r o p o r t i o n a l to the m u l t i d i m e n s i o n a l Gaussian probability density function
(
p(x,m,d)=exp
- ~
1,
(5)
where m, is a mean over the parents x , and d~ are the components o f a dispersion vector adjusted at each iteration. We have a (/2+2) strategy if in the next iteration the new parents are the # vectors with the smallest errors in the set comprising both parents and children. We have a (#, 2) strategy if the best # vectors are selected among the " c h i l d r e n " only ( 2 > #). The so-called Monte Carlo m e t h o d is an example of a ( 1 + 1 ) strategy. Fig. 1 graphically illustrates the (3 + 4) and the (2, 6) strategies. The other ingredient o f an evolution strategy is the rule by which the dispersions are modified; if the parents are close to the m i n i m u m and the range of search is much larger than the distance from it, only a small fraction of children will be better than the parents, and Idl should be reduced; when children are frequently better than parents the dispersion should be increased, to explore regions farther away from m. It follows that a decision about the dispersion at the next step should be taken on the basis o f a statistically significant n u m b e r of generations. 3
~=4
,L
£
g
p
X
X
X
- - >
(2)
Y
_
__
Y
Y
Y
(a) ~t 2
L=6 p
(I)
x
x
> Y
Y (.J) Y
YjLY)
(2)
X
X
> Y
Y
Y
Y
Y
Y
(b) Fig. 1. Graphic illustration of the ( 3 + 4) and the (2, 6 ) strategies.
Tuning elements of an evolution strategy are, typically, the length of the evolution histoD,, the fractional change of dispersion, the critical success rate o f children which triggers an increase of the search range. With increasing the dimension o f the x-domain we are likely to need to increase the evolution history,, but total c o m p u t a t i o n time may increase much less. The r a n d o m search explores simultaneously all x-dimensions, and the o p t i m i z e r overhead is almost independent of the n u m b e r of parameters. These stochastic methods have many attractive features, when c o m p a r e d with deterministic ones. Firstly, they are more effective in reaching a global minimizer, or at least a very stable local minimizer: secondly, the accuracy of solution is high because it does not d e p e n d upon the accuracy of the gradients: thirdly, the method can be applied even if the derivative of the objective function does not exist: and finally the treatment of constraint conditions is usually much simpler than in HODOMs. From the timeefficiency point of view, some results show that stochastic methods are comparable with deterministic ones, and sometimes they work even better [8,9].
3. T e s t i n g t h e o p t i m i z e r s
No optimizer warrants the convergence to the absolute m i n i m u m o f f ( x ) , nor every' problem can be treated with an optimizer. It has been argued [ 10] that d e t e r m i n a t i o n o f the E F G from a second-order quadrupolar spectrum is one of these problems, since substantially different EFGs may produce very similar powder patterns. An " e l e m e n t a r y " N M R MAS profile can be calculated at a given field B by assigning the strength of the q u a d r u p o l a r interaction VQ, the asymmetry p a r a m e t e r t/, and the width A of a broadening function (a Gaussian in our case). We have tested the optimizers by using as "experimental" functions spectra c o m p u t e d for known values of VQ, ~1, A. By trial and error we found the move limit adj u s t m e n t factor (0.95) which ensures convergence of our H O D O M to the "right" parameters, no matter how the initial point is chosen. In the process we identified several of the potential ambiguities denounced by Snyder et al. [ 10 ]. However, there is no doubt that with noiseless and undistorted data a ro-
P. Cofrancesco et aL /Quantitative interpretation of NMR spectra of disordered solids
871
Table 1 Data on optimization runs with a 1 + 1 evolution strategy q
F-calls
Errors of solution
18ql
8~,Q(%)
0.00 0.00
0.120 0.034
0.154 0.000
0.03 0.16
0.24 0.05
109 270
0.33 0.33
0.044 0.053
0.380 0.310
10.80 14.00
0.67 0.88
135 127
0.66 O.66
0.010 O.O21
0.460 1.230
5.04 1.36
0.46 1.15
400 147
1.00 1.00
0.082 0.053
2.380 1.530
3.72 1.67
0.35 0.22
167 269
bust optimizer will eventually approach the true solution. A typical solution has a rms< 0.1% and requires about 103 function calls. Data about optimization runs with a 1 + 1 evolution strategy are summarized in table 1. Each optimization has been repeated two times for spectra computed at q=0, 0.33, 0.66 and 1. The table evidences the order-of-magnitude fluctuations of the accuracy of the solution expressed as rms deviation, and over a 300% variation in the number of function calls. This fact had to be expected, since both starting point and evolution depend upon the rolling of a dice, and the random number generator is casually seeded before each run. Runs with about 100 function calls are typical, but accuracy in determination of the parameters (in particular q) is not always satisfactory. One may tune the parameters of the optimizer to reach better solutions in longer runs; alternatively, one may repeat several times the same run to estimate experimentally how accurate the solutions are, and what is the sensitivity of the error function near the minimum to the various parameters.
4. Two examples With a robust and reliable optimizer at our disposal, our attention should focus upon building an appropriate model of the experimental spectrum. In a polycrystalline powder with the resonant nuclei all in equivalent positions we may use the single-site (UQ, r/, A)-model used above. A should be strongly con-
8A (%)
rms (%)
strained with independent estimates of the instrumental broadening, and all causes of distortions discounted from the experimental data, or inserted in the model. To simulate the glassy counterpart of this solid we have to compromise between a realistic picture of disorder and the simplicity of the mathematical description. We may begin with the simplest of the models of disorder, and progressively refine it until a "reasonable" agreement with the experiment has been reached. An example of this procedure is illustrated for the case of the ~IB MAS spectrum of vitreous B203, chosen because it was expected to consist only of borons bonded to three bridging oxygens (B3). In model A we have assumed that disorder just results in an additional broadening, of the same nature of the instrumental one, and lifted the constraint upon A. Our HODOM gives a A about three times the instrumental one, but the rms error is high, 9.09%. Model B assumes that disorder is described by a normal distribution of q values, constrained between 0 and 1, centered at qo and with a width 2At/ at l / e P(qo). However, a substantially better agreement is reached (model C, fig. 2) by normally distributing the quadrupolar coupling constant PQ in a way statistically independent of the q-distribution:
P("Q)
L \ ~z,Q ] =]A
(6)
The broadening reaches the more reasonable value of 164 Hz, and the error reduces to 5.06%. However,
872
P. Cofrancesco et al. /Quantitative interpretation oJ'NMR spectra of disordered solids
tO moisture contamination through the reaction 2 B 3 + H20,--.2 ( B 4 - + H + ) ,
2,83
KHz
-2.88
Fig. 2. Model C. the difference plot Fe(Og) - F c ( o g ) (fig. 2) clearly indicates that the experimental spectrum contains small and narrow components, which we model with a Lorentzian-Gaussian profile 1
{l-g[
( ~ ) ]
2
--1
2
+0.69 g e x p [ - 0 . 6 9 ( ~o-)--- )e)o
]},
(7)
and four additional parameters: p, relative amplitude of the narrow component; ~Oo, position of the baricentre; g, relative weight of the Gaussian component; A', half height at half width of L + G functions. The parameters relative to the B3 sites remain essentially equal to those of model C, but the rms error reduces to 2.41% (model D). Position and width of the narrow component agree with the values characteristic of the component arising from fourcoordinated borons (B4). If this component is due
its amplitude ( p = 3%) allows to evaluate a reasonable water content of 0.7 wt%. Table 2 lists optimum values of parameters for models A, B, C and D. Fig. 3 shows an ~B MAS spectrum (open circles) obtained in the Li + glassy conductor Li20:B203: Li2SQ. The modelled one (solid curve) has been obtained with our H O D O M by assuming that four different types of structural units are present. (i) Symrnetric B3 units, with the same profile of the B3 units in B203. The only adjustable parameters of this profile are its weight p(B3S) and chemical shift. (ii) B4 6nits, with a small quadrupolar broadening, and a profile described by eq. (7). Beside the weight factor p(B4), the adjustable parameters are position (~o0) , width (A'), and fraction (g) of the Gaussian character of the profile. (iii) Asymmetric B3 units, with borons connected to one or two non-bridging oxygens. In an extensive N M R study of lithium-rich borate glasses, Yun and Bray [ l 1 ] have been able to determine average values Of UQand ~/for this unit; we have derived its EFG parameters by fitting our MAS spectra o f L i 2 0 : B 2 0 3 glasses. We obtained UQo= 1.38 MHz, and qo=0.43, values consistent with the values quoted by Yun and Bray [ 11 ]. In addition, we found a very small dispersion for vo ( A u o = 3 0 kHz), and a very large one for r/(At/= 1.96). As for the (i) component, the only adjustable parameters for the B3A profile are weight, p ( B 3 A ) , and chemical shift. (iv) Borate units bonded to sulphur, SB, with weight p(SB) = 1 - p ( B 3 S ) - p ( B 4 ) - p ( B 3 A ) . Their possible structures have been discussed in a previous
Table 2 Optimum values of the parameters for A, B, C and D models, Model
qo
Aq
UQo(MHz)
AUQ
A (kHz)
rms (%)
A B C
0.22 0 0
1.799 1.860
1.357 1.357 1.357
0.043
0.327 0.278 0.164
9.09 7.99 5.06
D
coo (kHz)
g
A' (kHz)
p(B4)
rms (%)
0.808
0.89
0.337
0.03
2.41
P. Cofrancesco et al. /Quantitative interpretation of NMR spectra of disordered solids
6.71
I(]-Iz
2.16
Fig. 3. ~B MAS spectrum (see text). Table 3 Optimal parameters p (B3S)=0.151 p (B3A)=0.145 p (B4) =0.416 090=4.908 kHz 4' =0.208 kHz g=0.659
p (SB) =0.288 t/o= 0.6 A~/= 2.89 UQo= 1.148 MHz AvQ= 1.52 zi=0.194 kHz
work [12]. Here, we estimate the EFG parameters of SB by letting the optimizer search for the best fit values of VOo, Az,o, qo and Aq, with essentially no constraints. With this model, a good agreement with experiment, and a low rms error of 3.54% is reached. Table 3 lists all data obtained from this optimization run, with the exclusion of the chemical shift parameters. We acknowledge that this model requires the checking of several assumptions (not discussed here) and implies the knowledge of complex profiles obtained from independent measurements; however, all this effort is justified by the rich set of structure-related parameters we have been able to extract from a single spectrum.
6. Conclusions A primary goal of the paper is that of conveying the idea that profile analysis of quadrupolar spectra is a daunting but rewarding task. With an ideal ex-
873
perimental setup we may plan MAS experiments at very high polarizing fields, which allow estimating number and amplitudes of the distinct components by making the chemical shift interaction predominant, and at relatively low magnetic fields, where we are more sensitive to the EGF distribution of the individual components. However, the second-order quadrupolar interaction is always associated with an overall displacement of the line, which we handle by simultaneously optimizing EFG and chemical shift parameters. Consistency of the latter ones over families of related glasses is a strong evidence in favour of the model chosen. When it can be proven that a family of glasses consists of a small number of structural units with characteristic EFG distributions, the analysis simplifies enormously, since we can adopt as spectral components the characteristic profiles of these units, and adjust their amplitudes to fit the experimental data. In general, however, a model results from linear combination of many primitive profiles, each representing the powder pattern corresponding to a choice of r/, u, and a. The computation time may become an issue even with the most powerful hardware currently available. We deal with it by selecting and tuning the optimization technique best suited to a given problem, and by writing a most efficient code for the error function. In N M R spectral modeling, the technique used to compute this function plays a very critical role. For this reason, we have evaluated several general methods for calculating, with prescribed resolution and accuracy, the "powder averages" leading to a primitive profile. The results of these studies will be presented elsewhere.
References [ 1 ] G. Engelhardt and D. Michel, High Resolution Solid State NMR of Silicates and Zeolites (Wiley, New York, 1987). [2 ] G.L. Turner, R.J. Kirkpatrick, S.H. Risbud and E. Oldfield, Ceram. Bull. 66 (1987) 656. [3] P. Mustarelli, M. Rudnicki, A. Savini, F. Savoldi and M. Villa, Mag. Res. Imag. 8 (1990) 101. [4] S. Kirkpatrick, C.D. Gelatt and M.P. Vecchi, Science 220 (1983) 671. [5] I.O. Bohachevsky, M.E. Johnson and M.L. Stein, Technometrics 28 (1986) 209. [6] L. Davis, Genetic Algorithms and Simulated Annealing (Morgan-Kauffmann, Los Altos, 1989 ).
874
P. Cofrancesco et al. / Quantitative interpretation qf NMR spectra qf disordered solids
[7] A. Gottvald, K. Preis and M. Friedrich, OptimiZation Methods for Computational Electromagnetics, IGTE Syrup., Graz (1990) p. 23. [8] A. Gottvald, K. Preis, C. Magele, O. Biro and A. Savini, Global Optimization Methods for Computational Electromagnetics, in: Proc. Compumag, Sorrento, Italy, Vol. 2 (1991) pp. 1077-1080.
[9] J.H. Kalivas, N. Roberts and J.M. Sutter, Anal. Chem. 61 (1989) 2024. [10] L.C. Snyder, G.E. Peterson and C.R. Kurkjian. J. Chem. Phys. 64 (1976) 1569. [ 11 ] Y.H. Yun and P.J. Bray, J. Non-Cryst. Solids 44 ( 1981 ) 227. [ 12 ] P. Mustarelli, S. Scotti, M. Villa and P.R. Gandhi, Solid Stale lonics 39 (1990) 217.