- Email: [email protected]
A New Method for Simulating Randomly Oriented Powder Spectra in Magnetic Resonance: TheSydneyOperaHouse(SOPHE) Method
JOURNAL OF MAGNETIC RESONANCE,
Series A 117, 1–8 (1995)
A New Method for Simulating Randomly Oriented Powder Spectra in Magnetic Resonance: The Sydney Opera House (SOPHE) Method DEMING WANG AND GRAEME R. HANSON Centre for Magnetic Resonance, University of Queensland, St. Lucia, Queensland, 4072, Australia Received October 25, 1994; revised March 29, 1995
A new method named the Sydney Opera House (SOPHE) method for computer reconstruction of randomly oriented powder spectra in magnetic resonance is presented. This method features a new way of partitioning the unit sphere in conjunction with global cubic spline interpolation and local linear interpolation. This allows randomly oriented powder spectra of high quality to be produced in significantly shorter times. This new method is especially efficient in simulating powder spectra of spin systems of large dimensions. It is a general method and is applicable to all types of randomly oriented powder spectra in magnetic resonance. q 1995 Academic Press, Inc.
The main battle to combat in computer reconstruction of randomly oriented powder spectra is anisotropy. If a system has a high degree of anisotropy, as often found in orthorhombic or lower symmetries, the randomly oriented powder spectrum must be reconstructed from a very large number of single-crystal spectra. Otherwise, the simulated spectrum shows considerable ripples, computer noise. For example, a rhombic Mn(II) complex with zero-field splittings comparable with the microwave frequency will yield a spectrum spanning a field range of Ç400 mT. Assuming the full widths at half maximum (FWHM) are on the order of 2 mT and that the tolerance of the resonant field position is not to exceed FWHM, a minimum of 200 steps in u will be required to cover 907 on an average base. By using the Igloo partition method (14), discussed later, some 20,000 orientations in the first octant must be sampled. Normally, a few diagonalizations, say an average of five, are required to obtain the resonant-field position for a particular transition for each orientation. For 30 allowed transitions, six from each hyperfine set, there will be a total of three million matrix diagonalizations required for the generation of the randomly oriented spectrum. Given that a single matrix diagonalization takes Ç0.1 seconds on a modern computer like a Sun SPARCstation 10/30 workstation, it would thus take 3.7 days to produce a single simulated spectrum. This is obviously impractical. One alternative to avoid a large number of single-crystal spectrum calculations is through interpolation, and this has been considered by a number of authors (6, 7, 10, 15). In these interpolation methods, the eigenvalues and eigenvectors are calculated exactly for only a small number of data points, while for the remaining data points, an interpolation method is used. This automatically reduces the computing time. A number of interpolation methods have been used in conjunction with different partition schemes (6, 7, 10, 15). Although these methods have been successfully applied to certain specific problems, there is limited generality to these methods. In this paper, we will present a new method for simulating
INTRODUCTION
Computer simulation has become a common method in magnetic resonance for extracting the spin Hamiltonian parameters from randomly oriented powder spectra or spectra of disordered spin systems (1, 2). However, this practice is still largely limited to relatively simple systems of low dimensions in spin space where analytical solutions exist or perturbation theory can be employed ( 2–5). A typical example is a spin S Å 12 interacting with a number of nuclear spins where the electronic Zeeman interaction dominates other interactions involved and the latter are treated as a second-order perturbation. For spin systems with S § 12, perturbation theory is generally not applicable because of the presence of other equally important interactions such as the fine structure. For a system in which a number of interactions with comparable magnitudes are present, the only answer to the problem seems to be matrix diagonalization; that is, all of the interactions are treated with equal importance and the energy matrix is diagonalized in the entire spin space (6–11). This approach can easily become a formidable task even in the modern era of computers as diagonalization of Hermitian matrices is a O(N 3 ) process, where N is the order of the matrix. For example, with a Sun SPARCstation 10/30 workstation, diagonalization of a 36 by 36 Hermitian matrix for both eigenvalues and eigenvectors with one of the fastest routines available, hsh (12) or lapack (13) routines, takes Ç0.1 seconds. 1
m3839$0679
10-10-95 22:18:52
magas
1064-1858/95 $12.00 Copyright q 1995 by Academic Press, Inc. All rights of reproduction in any form reserved.
AP: Mag Res
2
WANG AND HANSON
FIG. 1. (a) The Igloo method for partitioning an octant of a unit sphere (14). There are N values of u but the number of values for f is increased from 1 to N in steps of 1. Here N Å 10 (55 vertex points). (b) The triangular method used by Alderman et al. (15). In this method, the unit sphere (only the first octant shown here) is partitioned into identical equilateral triangles on the faces of an octahedron.
randomly oriented powder spectra. This method features a new way of partitioning a unit sphere into triangles of nearly equal solid angles, and a new interpolation scheme which includes a global cubic spline interpolation for data points on the edges of the triangles and linear interpolation for points inside the triangles. The efficiency and reliability of this new method are tested on a number of representative spin systems in EPR. PARTITIONING OF A UNIT SPHERE
Before our new partition method is introduced, a short review of the various partition methods which have been used in computer simulations is presented. The Apple-peel method is the simplest method to partition a sphere which involves the division of the polar angle u from 07 to 907 in N equal steps and the azimuthal angle f in M equal steps. Although this method is simple from the viewpoint of implementing an interpolation scheme (6, 7), it is wasteful as half of the data points lying between u Å 07 and 457 cover less than 30% of the total solid angle of the first octant. A modification of the apple-peel method has also been used by some authors (2). The Igloo method, developed by Belford et al. (14) partitions the unit sphere more evenly. An example of this method for N Å 10 is shown in Fig. 1a. There are N values of u but the number of values for f increases from 1 at u Å 07 to N at u Å 907. Thus, the total number of orientations required in the Igloo method is N(N / 1)/2 and is consequently computationally more economical than the apple-peel method (N 1 M). However, this method is unsuitable for the incorporation of interpolation algorithms. The spiral method was introduced by Mombourquette and Weil (10) and involves setting up a spiral circling about the polar axis (the z axis) of the unit sphere, starting at the x
m3839$0679
10-10-95 22:18:52
magas
axis and spiralling up to the pole. The central point of the spiral method is to keep the arc length between any consecutive points along the spiral constant and also keep the separation between the loops along any f constant. The disadvantage of the spiral method is the lack of an analytical solution for determining the next point from a given point, which must be determined from a minimization algorithm. Another method used in computer simulation of randomly oriented powder spectra, the triangular method, was introduced by Alderman, Solum, and Grant (15). It involves an equilateral triangular grid on the faces of an octahedron as shown in Fig. 1b. The plane triangles on the faces of the octahedron have the same area, but they do not all subtend the same solid angle. An exact calculation of the solid angles by individual triangles is complicated and thus only an approximation was given. Our new method, the Sydney Opera House (SOPHE) method, can be simply described by defining the vertex points on the unit sphere in exactly the same way as in the Igloo method but connecting the points as in the triangular method. Figure 2a shows an N Å 11 partition where N has the same meaning as in the Igloo method. There are (N 0 1) 2 triangles, the same number as found in the triangular method. Partitioning the unit sphere in this way resembles the roof of the famous Australian Sydney Opera House and consequently we have named our method SOPHE. The advantage in using SOPHE to partition the unit sphere lies not in the fact that the triangles in SOPHE subtend nearly the same solid angle but in its suitability for applying a new and powerful interpolation scheme which is the main subject of this paper. Before this new interpolation scheme is introduced, a closer examination of the grid employed in SOPHE is presented.
FIG. 2. The Sydney Opera House (SOPHE) method for partitioning a unit sphere. (a) The first octant revealing the grid defined by Eq. [1] in the text. (b) The angular grid produced by subpartitioning each triangle into ‘‘tiny’’ triangles.
AP: Mag Res
3
SIMULATING RANDOMLY ORIENTED POWDER SPECTRA
Similar proof for the down triangles is provided by the equation
THE GRID EMPLOYED IN SOPHE
It is not too difficult to see that the grid in Fig. 2a is made up of three sets of curves which are defined u Å C1
[1a]
uf Å C2
[1b]
uf Å
pu 0 C2 . 2
[1c]
The most obvious set of curves, set 1, are those which have constant u (Eq. [1a]). C1 takes on (N 0 1) equally spaced values, C1 Å iDu, where i Å 1, 2, . . ., N 0 1 and Du Å p /[2(N 0 1)]. Less obvious are the other two sets of curves. Curves governed by Eq. [1b] are referred to as set 2, and C2 is also a constant equal to (i 0 1) Dup /2 with i taking the same set of values as for C1 . When i Å 1, it gives a curve for f Å 07 in the xz plane, the only curve which needs to be considered for the simulation of axially symmetric, randomly oriented, powder spectra. The remaining set, set 3, is defined by Eq. [1c]. Clearly, the grid only depends on one parameter, N, which determines Du. The mathematical realization of the grid in SOPHE made it possible to derive the analytical expressions for calculating the solid angles subtended by individual triangles. In order to facilitate this derivation and discussions later on in the interpolation scheme, we introduce a few definitions. We define a band as consisting of all triangles between u1 Å (i 0 1) Du and u2 Å iDu. For a given partition number N, there are N 0 1 bands. For a given band, except for the first one, which has only one triangle, there are two types of triangles, called up ( n ) and down ( , ) triangles. The number of up triangles is equal to i and the number of down triangles is one less, i 0 1. We now show that, in a given band, all up triangles subtend the same solid angle as do all the down triangles. With the help of Eq. [1b] and Eq. [1c], the analytical expression for the solid angle for each of the up triangles in the band defined by u1 Å (i 0 1) Du and u2 Å iDu is given as d vD Å 0
*
u2
d cos u
u1
*
( p / 2 ) 0 ( p / 2 ) Du( i0p ) / u
df
( p / 2 ) Du( p 01 ) / u
Å
p (cos u1 0 cos u2 ) 2
0
p (i 0 1)Du 2
*
u2
u1
sin u du, u
[2]
where p is the up triangle index taking values from 1 to i. As shown, d vD is independent of the index p, thus proving that all of the up triangles subtend the same solid angle.
m3839$0679
10-10-95 22:18:52
magas
d vÇ Å 0
*
u2
d cos u
/
df
( p / 2)0 (p / 2 ) Du( i0p ) / u
u1
Å0
*
( p / 2 ) Du p / u
p (cos u1 0 cos u2 ) 2 p iDu 2
*
u2
u1
sin u du. u
[3]
Now p takes i 0 1 values from 1 to i 0 1 as there are only i 0 1 down triangles in the band. Evidently, d v, is also independent of the index p. The total solid angle subtended by the given band is equal to id vD / (i 0 1)d vÇ Å (p /2)(cos u1 0 cos u2 ), as expected. Both d vD and d v, involve the evaluation of the integral (sin u / u ) du. This integral can be expanded into a fast converging series so that it can be easily evaluated numerically. However, for typical partition numbers (N) used in the simulation of randomly oriented powder spectra, the difference between d vD and d v, is very small. For example, when N Å 19 ( Du Å 57 ) dv, 0 d vD/d v, is 1.76 1 10 02 for the band defined by u1 Å 857 and u2 Å 907, the band in which the largest difference is found. This difference diminishes to 1.27 1 10 03 for the band defined by u1 Å 57 and u2 Å 107, the band of the smallest difference. With a larger N, the difference becomes smaller, which is as expected. Taking N Å 91 as another example, ( Du Å 17 ), the respective values are 3.7 1 10 03 and 5.08 1 10 05 . Each triangle in Fig. 2a can be easily subpartitioned into smaller triangles, referred to as tiny triangles. In Fig. 2b, a selected triangle is further partitioned into 81 tiny triangles with a subpartion number M Å 10. The grid formed in such a subpartition can still be described by Eq. [1]. In this particular case, C1 is stepped in a smaller step of p /[2(N 0 1)(M 0 1)] from C1 Å 457 to C1 Å 547, the two corresponding curves which bound the triangle. A similar process is applied to curves in sets 2 and 3. The subpartition is used in conjunction with a linear interpolation scheme. Although we can calculate the solid angles corresponding to each of these tiny triangles by using Eq. [2] and Eq. [3], such exact calculations are hardly justifiable in view of the approximations involved in the linear interpolation. For example, for the N Å 19 partition with a subpartition number M Å 51 (in other words, each triangle is subpartitioned into 2500 tiny triangles), the biggest difference among the solid angles subtended by the tiny triangles in a given triangle was found to be less than 5%. Therefore, in our program, we assumed all the tiny triangles in a given triangle subtend the same solid angle, thus saving many calculations which would be required if Eq. [2] and Eq. [3] were used.
AP: Mag Res
4
WANG AND HANSON
A NEW INTERPOLATION SCHEME
Most interpolation schemes adopted so far are two-dimensional (6, 7, 15) with u and f as variables. One-dimensional interpolation has been considered by Mombourquette and Weil (10) but they only considered the simplest interpolation method, linear interpolation. In most cases, linear interpolation cannot satisfactorily reconstruct the experimental spectra without performing exact calculations of the positions and intensity at a very large number of orientations. When higher-order interpolation methods are considered, first derivatives with respect to u and f are often required. The task of calculating these first derivatives can easily become as formidable as obtaining the position/intensity information. An important feature in our new interpolation scheme is that we only use one-dimensional interpolation methods in which the demand for the first-derivative calculations is dramatically reduced. Specifically, we use cubic spline interpolation for data points along the edges of the triangles and linear interpolation for those inside the triangles from the data obtained from the cubic spline interpolations. Global Cubic Spline Interpolation Cubic spline interpolation is featured by its smoothness in the first derivative and continuity in the second derivative, both within an interval and at its boundaries. Therefore, it should be a good interpolation method for the data points along the edges of the triangles. Strictly speaking, among the three sets of curves in SOPHE, only set 1 (Eq. [1a]) involves one variable, f. The other two sets of curves involve both u and f. However, by using Eq. [1b] and Eq. [1c], sets 2 and 3 can also be treated as involving only one variable. In both sets, we keep u as the variable and treat f as a parameter and thus one-dimensional interpolation methods can be used for all three sets. The cubic spline interpolation process is divided into two steps (16). For a given set of data, xi and yi Å y(xi ), i Å 1, . . ., n, and two end conditions, the first step is to calculate the second derivatives d 2 y/dx 2 (or the first derivatives, dy/ dx) at the points xi (i Å 1, 2, . . ., n). Different types of end conditions can be used (16) and here we adopt the firstderivative end condition, which is to specify the two first derivatives at the two endpoints x1 and xn . The second derivatives are then obtained with a C function called spline in ‘‘Numerical Recipes in C’’ (17). For a given value x, the second step is to find the interval within which the specified x is located. In the algorithm we developed, this location search has been avoided, thus saving some computing time. In the actual calculation, we store the second derivatives. There are 6(N 0 1) first derivatives, 2(N 0 1) for each set of interpolation. This number is significantly reduced from what is required in the two-dimensional schemes used by
m3839$0679
10-10-95 22:18:52
magas
previous authors (6, 7). For example, in the interpolation scheme used by Gribnau, van Tits, and Reijerse (7), a total of 3N 2 first derivatives were required. For a partition number of N Å 31, it requires 16 times as many first-derivative calculations as used in our scheme. For each set, there are N 0 1 cubic spline interpolations with the number of data points from 2 to N in steps of 1. We have tested the quality of the interpolated data from the cubic spline interpolation by comparing them with those obtained through exact calculations on a number of spin systems and the results were found to be excellent. Here, we examine a specific example of an S Å 32 spin system. The spin Hamiltonian used has the following form K Å g bBrS / D [S 2z 0 13 S(S / 1)] / E(S 2x 0 S 2y ),
[4]
where the terms have their usual meanings (2). The comparison was made on a highly anisotropic system in which E/D Å 13, the so-called fully rhombic symmetry. We used the data for a Cr(CN) 30 center in a KCl lattice (18). For 6 this center, D Å 2.73 GHz, g à 1.990, and E/D Å 13. The microwave frequency was set to 9.5 GHz. We made the comparison on the three allowed transitions, 1 } 2, 2 } 3, 3 } 4 with the energy levels numbered in order of increasing energy. In the exact calculation, we used a bisection method in locating the resonant positions with an accuracy of {0.01 mT. It should be pointed out that the bisection method is one of the least efficient methods in locating resonant fields and is used simply for the purpose of demonstrating the SOPHE algorithm. The resonant field position of the allowed transitions varies between 160 and 536 mT. Table 1 presents a summary of the comparison. It also includes the maximum and minimum resonant fields for each transition. The data points selected for comparison were obtained by stepping u by 0.57 for curve sets 2 and 3, while points of separation with approximately equal arc lengths of 0.57 were selected for curve set 1. The overall agreement in Table 1 is excellent. Naturally, when N gets larger, the grid becomes finer and the agreement is better. With N Å 31, no data points with differences between the calculated and interpolated larger than 0.02 mT were found for transition 3 } 4 and only 12 such points were detected for transitions 1 } 2 and 2 } 3, which represents only 0.16% of the total number of points compared. The average deviation is very small, a few microtesla for all three transitions. A similar comparison can be made for the intensity but normally the transition probability varies much less than the resonant field positions and so this was not included in Table 1. Another important feature worth mentioning explicitly about Table 1 is that normally one would expect a better agreement for transition 2 } 3 as its resonant field varies
AP: Mag Res
5
SIMULATING RANDOMLY ORIENTED POWDER SPECTRA
TABLE 1 Comparison of the Resonant Fields Interpolated through the Global Cubic Spline Interpolation (Binterpol) with Those Calculated by Full Matrix Diagonalization (Bexact) N
a
dub (degrees)
Npc
3}4
Bmaxg Å 536.0039 (y axis)
19 25 31
0.01 0.01 0.01
1}2
Bmax Å 536.0039 (z axis)
19 25 31
0.01 0.01 0.01
2}3
Bmax Å 330.5821 (y or z axis)
19 25 31
0.01 0.01 0.01
4842 6488 7480
4842 6488 7480
4842 6488 7480
K1d
K2e
Average deviation f (mT)
Local Linear Interpolation
Bmin Å 160.7016 (z axis) 0 0 0
3 0 0
0.90 0.31 0.16
Bmin Å 160.7016 (y axis) 19 9 0
71 26 12
1.97 0.69 0.33
Bmin Å 299.7798 (x axis) 19 7 0
70 25 12
1.68 0.59 0.29
a
N is the partition number. du is the step used for evaluating the first derivatives. c Np is the total number of data points compared. d K1 is the number of data points for which DB § 50 mT (DB Å ÉBinterpol 0 BexactÉ). e K2 is the number of data points for which DB § 20 mT. f Average deviation is defined as the average of DB for all the points compared. g Bmax (mT) and Bmin (mT) are the maximum and minimum resonant field positions when the magnetic field is applied along the axes shown in brackets. b
only in a small fraction of the field range spanned by the other two transitions, but the opposite is shown when comparing transitions 2 } 3 and 1 } 2. The reason for this is that the quality of the interpolated data depends on the accuracy of the first derivatives. As said, the first derivatives in the full diagonalization method must be obtained numerically, and to evaluate all the first derivatives to a good precision is not an easy task. This is because first derivatives normally take values over a wide range. In our test case, they vary from a few hundred microtesla/degree to over 10 mT/degree, and as they are not known a priori, it is impossible to choose a good value for du or df for each firstderivative evaluation. In our program, we set the same value for du for all first-derivative calculations. For curve sets 2 and 3, an appropriate amount of df is necessary, which is obtained through Eq. [1b] and Eq. [1c]. For curve set 1, we choose values for df by using df Å du/sin u. We have
m3839$0679
10-10-95 22:18:52
tested the quality for the interpolated data as a function of du for all three transitions and the results for transition 1 } 2 are shown in Table 2. In general, any values between 0.0057 and 0.17 are acceptable. However, we found 0.017 is an optimal value for du for our test case with the resonant field positions determined to an accuracy of {0.01 mT.
magas
After the global interpolation, the position/intensity/linewidth information can be viewed as a ‘‘skeleton’’ based on the designed grid in Fig. 2a. Further subpartitioning can be made for each triangle, and for any data point inside the triangle, we employ a linear interpolation method. This process can be viewed as a ‘‘tile-filling’’ process. One advantage of the linear interpolation method is that it does not involve any derivatives, so that no further exact calculations are required. The justification of the linear interpolation relies on the anisotropy of the system under consideration and such a justification should always be achievable by varying the grid size. A number of linear interpolation schemes can be employed, which are schematically shown in Fig. 3. The singleset scheme shown in Fig. 3a suits most situations where the resonant field position hardly changes along one of the edges. In this case, we linearly interpolate the data points with the set which includes the edge along which the resonant field position shows little change. For other situations, the combined (Fig. 3b) or the averaged scheme (Fig. 3c) can be used. Obviously, the averaged scheme requires more computing time. However, we noticed that the averaged scheme in general produced higher-quality (less ripple) simulations
TABLE 2 Comparison of the Resonant Fields Interpolated through the Global Cubic Spline Interpolation Method (Binterpol) with Those Calculated by the Full-Matrix-Diagonalization Method (Bexact) for Transition 1}2 with Different Values of dua N
du (7)
Np
K1
K2
Average deviation (mT)
31 31 31 31 31 31 31
0.005 0.01 0.05 0.1 0.5 1.0 3.0
41,290 41,290 41,290 41,290 41,290 41,290 41,290
0 0 0 7 496 2632 6902
57 55 66 116 2091 4873 10058
0.33 0.31 0.51 0.82 3.50 53.3 98.0
a The symbols have the same definitions as in Table 1. Here, a smaller step of 0.17 (0.57 was used for Table 1) was used in selecting the points for comparison in order to show there was absolutely no oscillations in the interpolation. Consequently, the total number of data points increases from 7480 to 41,290 for N Å 31. The vertex points were excluded in all comparisons.
AP: Mag Res
6
WANG AND HANSON
FIG. 3. Schematic representations of the possible linear interpolation schemes. (a) Single-set scheme. (b) Combined scheme. (c) Averaged scheme.
than the combined scheme. Thus, in our program, the averaged scheme is used. In the averaged scheme, each point inside the triangle is interpolated three times from three different sets and an averaged value is taken for each point. RECONSTRUCTION OF A RANDOMLY ORIENTED POWDER SPECTRUM
In this section, we present some details of the computer implementation of SOPHE. Suppose the position/intensity/ linewidth information for a selected transition has been obtained at the vertex points for a given partition number N through a certain method. We can then start to calculate the first derivatives required for the cubic spline interpolation. If only the resonant position is globally interpolated, there are 6(N 0 1) first derivatives, 2(N 0 1) for each set of global interpolations. If both resonant position and intensity are globally interpolated, the number of first derivatives is simply doubled. Second derivatives are obtained with the function spline (17). These second derivatives, together with the position/intensity/linewidth information at the vertex points, are passed on to a function which handles interpolation and convolution. The interpolation process is made triangle by triangle and band through band. It starts with band i Å 1 at the top of the grid and progresses toward the xy plane. For a given band, the data elements corresponding to the vertex points can be arranged as a one-dimensional array
m3839$0679
10-10-95 22:18:52
magas
so that to move from one triangle to the next is simply to shift the index by 1. Figure 4 shows such an indexing scheme for a band with five triangles. The three vertex points which define a triangle are j, j / 1, j / 2 with j shifting from 1 to 5. Therefore, the whole interpolation/convolution process involves only repeating the same routine which handles a single triangle and that is how we implemented the SOPHE method. In the actual calculation, we first obtain a profile of wkl d v against the resonant position for a selected triangle; wkl represents the intensity for the transition labeled by kl. This is done by going through all the tiny triangles into which the given triangle is partitioned. Then the profile obtained is convoluted with a given lineshape with the linewidth averaged from those at the vertex points before adding its contribution to the final spectrum. This process is repeated for every transition. Since EPR is a field-swept technique, corrections for the anisotropic transition probability and linewidth must be dealt with in a proper way in computer simulations. This is especially important for high-spin systems. When perturbation theory is used, these corrections can be taken into account point by point. This normally involves calculation of the spectrum in frequency space and then transforming back into field space and requires very little additional computational time (2). However, in the full-matrix-diagonalization method, such point-by-point corrections become impractical. We have, instead, used an approach introduced by van Veen (6) in which the corrections for the anisotropic transition probability and linewidth are applied at the triangle level, i.e., the vertex points where exact calculations are performed. An S Å 52 spin system with g Å 1.990, D Å 1.0 GHz and E/D Å 13 is used to demonstrate the efficiency of this algorithm based on the SOPHE method. Here, again we assume full rhombic symmetry with a large zero-field splitting, 7.06 GHz, which is comparable with the microwave frequency used, 9.0 GHz. A partition number N Å 19 was used, which means exact calculations were performed at 190 vertex points. In Fig. 5a, the interpolation was made in steps of 17, corresponding to 3905 interpolated points per transition for each of the five allowed transitions. As shown, the reconstructed spectrum is very ‘‘noisy.’’ In Fig. 5b, the interpolation was made in steps of 0.17, resulting in 405,260 data
FIG. 4. The indexing scheme used in our computer implementation of the SOPHE method.
AP: Mag Res
SIMULATING RANDOMLY ORIENTED POWDER SPECTRA
7
points per transition and the computer ‘‘noise’’ has been removed nearly completely. Figure 5b consumed 120 s more, 24 s/transition, in CPU time on our Sun SPARCstation 10/ 30 workstation. This means that at a marginal cost of 24 s, more than 400,000 data points have been generated. In fact, the large proportion of the consumed CPU time was spent in the sorting routine and the data-generating speed in the interpolation routine was about 50,000 data points per second, a remarkable achievement for the SOPHE method. Further improvement in spectrum quality can be made by using a finer grid (a larger value for N). In Fig. 5c, a partition number of N Å 31 ( Du Å 37 ) was used with a subpartition of 0.17. DISCUSSION AND CONCLUSION
5
FIG. 5. Simulated spectra for a model S Å 2 spin system (g Å 1.990, 1 D Å 1.0 GHz, and E/D Å 3). (a) Parameters used in the simulation were partition number, N Å 19 (190 vertex points); subpartition number, 6 (3905 interpolated points per transition); field axis resolution, 4096 points; field sweep-width, 350 mT; center field, 325 mT; microwave frequency, 9.0 GHz. A Gaussian lineshape with a half-width at half-maximum of 30 MHz was used. The linewidth was assumed to be isotropic. (b) Parameters used in the simulation were the same as for (a) except for a larger subpartition number of 51 (405,260 interpolated points per transition). (c) Parameters used in the simulation were the same as for (b) except with a finer grid, N Å 31 (496 vertex points).
m3839$0679
10-10-95 22:18:52
magas
The most significant features of this new algorithm based on the SOPHE method are its efficiency in generating a very large number of data points from those at vertexes, which are normally obtained through a time-consuming process, and its reliability. Data can be interpolated extremely rapidly, tens of thousands of data points per second on modern computers. For example, on our Sun SPARCstation 10/30, the data-generating speed in the interpolation process is Ç50,000 points/s. Thus, the SOPHE method has reduced the computing task mainly to calculating the position/intensity information at the vertex points, computing the first derivatives, sorting, and convoluting. This is a remarkable improvement in computer reconstruction of randomly oriented powder spectra. In particular, the SOPHE method should significantly ease the burden of simulating randomly oriented powder spectra of spin systems of large dimensions. The quality of the interpolated data has been demonstrated to be very high. The interpolated values on the edges can be expected to be very close to the true values if the first derivatives are computed with a good precision as typified in Table 1. This is important as the globally interpolated data are to form the base on which local linear interpolations are performed. Although we have only discussed the application of the SOPHE method to orthorhombic spin systems, it can be easily extended to monoclinic (C2 , C2h , and Cs ) and triclinic (C1 ) symmetries where two and four octants are required respectively. In fact, the SOPHE method introduced in the previous section is a special case of the more general situation where a section of the unit sphere is defined as u1 Å 07 and u2 Å 907, and f1 and f2 are two arbitrary angles with f2 ú f1 . Equation [1] can still be used providing f is replaced by f 0 f1 and p /2 is replaced by f2 0 f1 in both Eq. [1c] and the constant C2 . Therefore, for monoclinic symmetries, the SOPHE method can be applied by choosing f1 Å 07 and f2 Å 1807 if the z axis is defined as the symmetry axis (C2 and C2h ) or as the axis perpendicular to the reflection
AP: Mag Res
8
WANG AND HANSON
plane (Cs ). In the case of triclinic symmetry, half a sphere must be considered, which means f1 Å 07 and f2 Å 3607. In this case, set 1 curves are closed curves and periodic end conditions must be used. Alternatively, in monoclinic and triclinic symmetries, each individual octant can be treated separately in the same manner as for orthorhombic symmetry and their contributions are then added to the final spectrum. Subpartitioning of each triangle can also be easily made on an individual basis; that is, how fine each individual triangle is partitioned into tiny triangles is determined by the difference in the resonant field positions at the three vertex points. In other words, it can be determined by the ‘‘local’’ anisotropy for a given transition. Thus, the SOPHE method has a great flexibility in a sense that it handles the anisotropy step by step through the space, and, of course, transition by transition. The SOPHE method is a general method for simulating randomly oriented magnetic resonance spectra. It can be used for any computer reconstruction where the spectral function must be sampled at a very fine grid of the unit sphere, in other words, at a very large number of data points, and sampling at each point takes a considerable amount of computing time. For example, the method should find its use in simulating randomly oriented powder spectra in solid state NMR, electron nuclear double resonance (ENDOR), electron spin-echo envelope modulation (ESEEM), and nuclear quadrupole resonance in solids. Specific details concerning the availability of the SOPHE software can be obtained by writing to the authors.
REFERENCES 1. P. C. Taylor, J. F. Baugher, and H. M. Kriz, Chem. Rev. 75, 203 (1974), and references therein. 2. J. R. Pilbrow, ‘‘Transition Ion Electron Paramagnetic Resonance,’’ Clarendon Press, Oxford, 1990, and references therein. 3. G. M. Muha, J. Magn. Reson. 49, 431 (1982); 53, 85 (1983). 4. Y. Siderer and Z. Luz, J. Magn. Reson. 37, 449 (1980). 5. L. Andreozzi, M. Giordano, and D. Leporini, J. Magn. Reson. A 104, 166 (1993). 6. G. van Veen, J. Magn. Reson. 38, 91 (1978). 7. M. C. M. Gribnau, J. L. C. van Tits, and E. J. Reijerse, J. Magn. Reson. 90, 474 (1990). 8. G. G. Belford, R. L. Belford, and J. F. Burkhalter, J. Magn. Reson. 11, 251 (1973). 9. A. Kretter and J. Huttermann, J. Magn. Reson. 93, 12 (1991). 10. M. J. Mombourquette and J. A. Weil, J. Magn. Reson. 99, 37 (1992). 11. L. Cugunov, A. Mednis, and J. Kliava, J. Magn. Reson. A 106, 153 (1994). 12. M. Krusche, M. Neu, and H. Ziegler, J. Magn. Reson. 95, 368 (1991). 13. E. Anderson, Z. Bai, C. Bischof, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, S. Ostrouchov, and D. Sorensen, ‘‘LAPACK Users’ Guide,’’ Society for Industrial and Applied Mathematics, Philadelphia, 1992. 14. (a) M. J. Nilges, Ph.D. thesis, University of Illinois, Urbana, Illinois, 1979; (b) R. L. Belford and M. J. Nilges, Computer Simulation of Powder Spectra, EPR Symposium 21st Rocky Mountain Conference, Denver, Colorado, 1979; (c) A. M. Maurice, Ph.D. thesis, University of Illinois, Urbana, Illinois. 15. D. W. Alderman, M. S. Solum, and D. M. Grant, J. Chem. Phys. 84, 3717 (1986).
ACKNOWLEDGMENTS
16. Su Bu-qing and Liu Ding-yuan, ‘‘Computational Geometry— Curves and Surface Modeling,’’ Academic Press, San Diego, 1989.
We thank Doug Maclachlan and Kevin Gates for providing us some fast diagonalization routines and useful references, and Simon Teed for drawing some of the figures. The Australian Research Council is also thanked for awarding Research Grant A29331095 to G.R.H. which made this work possible.
17. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, ‘‘Numerical Recipes in C,’’ p. 113, Cambridge Univ. Press, Cambridge, 1992.
m3839$0679
10-10-95 22:18:52
magas
18. D. M. Wang, D. H. Hutton, and J. R. Pilbrow, J. Phys. C, 19, 789 (1986).
AP: Mag Res
Series A 117, 1–8 (1995)
A New Method for Simulating Randomly Oriented Powder Spectra in Magnetic Resonance: The Sydney Opera House (SOPHE) Method DEMING WANG AND GRAEME R. HANSON Centre for Magnetic Resonance, University of Queensland, St. Lucia, Queensland, 4072, Australia Received October 25, 1994; revised March 29, 1995
A new method named the Sydney Opera House (SOPHE) method for computer reconstruction of randomly oriented powder spectra in magnetic resonance is presented. This method features a new way of partitioning the unit sphere in conjunction with global cubic spline interpolation and local linear interpolation. This allows randomly oriented powder spectra of high quality to be produced in significantly shorter times. This new method is especially efficient in simulating powder spectra of spin systems of large dimensions. It is a general method and is applicable to all types of randomly oriented powder spectra in magnetic resonance. q 1995 Academic Press, Inc.
The main battle to combat in computer reconstruction of randomly oriented powder spectra is anisotropy. If a system has a high degree of anisotropy, as often found in orthorhombic or lower symmetries, the randomly oriented powder spectrum must be reconstructed from a very large number of single-crystal spectra. Otherwise, the simulated spectrum shows considerable ripples, computer noise. For example, a rhombic Mn(II) complex with zero-field splittings comparable with the microwave frequency will yield a spectrum spanning a field range of Ç400 mT. Assuming the full widths at half maximum (FWHM) are on the order of 2 mT and that the tolerance of the resonant field position is not to exceed FWHM, a minimum of 200 steps in u will be required to cover 907 on an average base. By using the Igloo partition method (14), discussed later, some 20,000 orientations in the first octant must be sampled. Normally, a few diagonalizations, say an average of five, are required to obtain the resonant-field position for a particular transition for each orientation. For 30 allowed transitions, six from each hyperfine set, there will be a total of three million matrix diagonalizations required for the generation of the randomly oriented spectrum. Given that a single matrix diagonalization takes Ç0.1 seconds on a modern computer like a Sun SPARCstation 10/30 workstation, it would thus take 3.7 days to produce a single simulated spectrum. This is obviously impractical. One alternative to avoid a large number of single-crystal spectrum calculations is through interpolation, and this has been considered by a number of authors (6, 7, 10, 15). In these interpolation methods, the eigenvalues and eigenvectors are calculated exactly for only a small number of data points, while for the remaining data points, an interpolation method is used. This automatically reduces the computing time. A number of interpolation methods have been used in conjunction with different partition schemes (6, 7, 10, 15). Although these methods have been successfully applied to certain specific problems, there is limited generality to these methods. In this paper, we will present a new method for simulating
INTRODUCTION
Computer simulation has become a common method in magnetic resonance for extracting the spin Hamiltonian parameters from randomly oriented powder spectra or spectra of disordered spin systems (1, 2). However, this practice is still largely limited to relatively simple systems of low dimensions in spin space where analytical solutions exist or perturbation theory can be employed ( 2–5). A typical example is a spin S Å 12 interacting with a number of nuclear spins where the electronic Zeeman interaction dominates other interactions involved and the latter are treated as a second-order perturbation. For spin systems with S § 12, perturbation theory is generally not applicable because of the presence of other equally important interactions such as the fine structure. For a system in which a number of interactions with comparable magnitudes are present, the only answer to the problem seems to be matrix diagonalization; that is, all of the interactions are treated with equal importance and the energy matrix is diagonalized in the entire spin space (6–11). This approach can easily become a formidable task even in the modern era of computers as diagonalization of Hermitian matrices is a O(N 3 ) process, where N is the order of the matrix. For example, with a Sun SPARCstation 10/30 workstation, diagonalization of a 36 by 36 Hermitian matrix for both eigenvalues and eigenvectors with one of the fastest routines available, hsh (12) or lapack (13) routines, takes Ç0.1 seconds. 1
m3839$0679
10-10-95 22:18:52
magas
1064-1858/95 $12.00 Copyright q 1995 by Academic Press, Inc. All rights of reproduction in any form reserved.
AP: Mag Res
2
WANG AND HANSON
FIG. 1. (a) The Igloo method for partitioning an octant of a unit sphere (14). There are N values of u but the number of values for f is increased from 1 to N in steps of 1. Here N Å 10 (55 vertex points). (b) The triangular method used by Alderman et al. (15). In this method, the unit sphere (only the first octant shown here) is partitioned into identical equilateral triangles on the faces of an octahedron.
randomly oriented powder spectra. This method features a new way of partitioning a unit sphere into triangles of nearly equal solid angles, and a new interpolation scheme which includes a global cubic spline interpolation for data points on the edges of the triangles and linear interpolation for points inside the triangles. The efficiency and reliability of this new method are tested on a number of representative spin systems in EPR. PARTITIONING OF A UNIT SPHERE
Before our new partition method is introduced, a short review of the various partition methods which have been used in computer simulations is presented. The Apple-peel method is the simplest method to partition a sphere which involves the division of the polar angle u from 07 to 907 in N equal steps and the azimuthal angle f in M equal steps. Although this method is simple from the viewpoint of implementing an interpolation scheme (6, 7), it is wasteful as half of the data points lying between u Å 07 and 457 cover less than 30% of the total solid angle of the first octant. A modification of the apple-peel method has also been used by some authors (2). The Igloo method, developed by Belford et al. (14) partitions the unit sphere more evenly. An example of this method for N Å 10 is shown in Fig. 1a. There are N values of u but the number of values for f increases from 1 at u Å 07 to N at u Å 907. Thus, the total number of orientations required in the Igloo method is N(N / 1)/2 and is consequently computationally more economical than the apple-peel method (N 1 M). However, this method is unsuitable for the incorporation of interpolation algorithms. The spiral method was introduced by Mombourquette and Weil (10) and involves setting up a spiral circling about the polar axis (the z axis) of the unit sphere, starting at the x
m3839$0679
10-10-95 22:18:52
magas
axis and spiralling up to the pole. The central point of the spiral method is to keep the arc length between any consecutive points along the spiral constant and also keep the separation between the loops along any f constant. The disadvantage of the spiral method is the lack of an analytical solution for determining the next point from a given point, which must be determined from a minimization algorithm. Another method used in computer simulation of randomly oriented powder spectra, the triangular method, was introduced by Alderman, Solum, and Grant (15). It involves an equilateral triangular grid on the faces of an octahedron as shown in Fig. 1b. The plane triangles on the faces of the octahedron have the same area, but they do not all subtend the same solid angle. An exact calculation of the solid angles by individual triangles is complicated and thus only an approximation was given. Our new method, the Sydney Opera House (SOPHE) method, can be simply described by defining the vertex points on the unit sphere in exactly the same way as in the Igloo method but connecting the points as in the triangular method. Figure 2a shows an N Å 11 partition where N has the same meaning as in the Igloo method. There are (N 0 1) 2 triangles, the same number as found in the triangular method. Partitioning the unit sphere in this way resembles the roof of the famous Australian Sydney Opera House and consequently we have named our method SOPHE. The advantage in using SOPHE to partition the unit sphere lies not in the fact that the triangles in SOPHE subtend nearly the same solid angle but in its suitability for applying a new and powerful interpolation scheme which is the main subject of this paper. Before this new interpolation scheme is introduced, a closer examination of the grid employed in SOPHE is presented.
FIG. 2. The Sydney Opera House (SOPHE) method for partitioning a unit sphere. (a) The first octant revealing the grid defined by Eq. [1] in the text. (b) The angular grid produced by subpartitioning each triangle into ‘‘tiny’’ triangles.
AP: Mag Res
3
SIMULATING RANDOMLY ORIENTED POWDER SPECTRA
Similar proof for the down triangles is provided by the equation
THE GRID EMPLOYED IN SOPHE
It is not too difficult to see that the grid in Fig. 2a is made up of three sets of curves which are defined u Å C1
[1a]
uf Å C2
[1b]
uf Å
pu 0 C2 . 2
[1c]
The most obvious set of curves, set 1, are those which have constant u (Eq. [1a]). C1 takes on (N 0 1) equally spaced values, C1 Å iDu, where i Å 1, 2, . . ., N 0 1 and Du Å p /[2(N 0 1)]. Less obvious are the other two sets of curves. Curves governed by Eq. [1b] are referred to as set 2, and C2 is also a constant equal to (i 0 1) Dup /2 with i taking the same set of values as for C1 . When i Å 1, it gives a curve for f Å 07 in the xz plane, the only curve which needs to be considered for the simulation of axially symmetric, randomly oriented, powder spectra. The remaining set, set 3, is defined by Eq. [1c]. Clearly, the grid only depends on one parameter, N, which determines Du. The mathematical realization of the grid in SOPHE made it possible to derive the analytical expressions for calculating the solid angles subtended by individual triangles. In order to facilitate this derivation and discussions later on in the interpolation scheme, we introduce a few definitions. We define a band as consisting of all triangles between u1 Å (i 0 1) Du and u2 Å iDu. For a given partition number N, there are N 0 1 bands. For a given band, except for the first one, which has only one triangle, there are two types of triangles, called up ( n ) and down ( , ) triangles. The number of up triangles is equal to i and the number of down triangles is one less, i 0 1. We now show that, in a given band, all up triangles subtend the same solid angle as do all the down triangles. With the help of Eq. [1b] and Eq. [1c], the analytical expression for the solid angle for each of the up triangles in the band defined by u1 Å (i 0 1) Du and u2 Å iDu is given as d vD Å 0
*
u2
d cos u
u1
*
( p / 2 ) 0 ( p / 2 ) Du( i0p ) / u
df
( p / 2 ) Du( p 01 ) / u
Å
p (cos u1 0 cos u2 ) 2
0
p (i 0 1)Du 2
*
u2
u1
sin u du, u
[2]
where p is the up triangle index taking values from 1 to i. As shown, d vD is independent of the index p, thus proving that all of the up triangles subtend the same solid angle.
m3839$0679
10-10-95 22:18:52
magas
d vÇ Å 0
*
u2
d cos u
/
df
( p / 2)0 (p / 2 ) Du( i0p ) / u
u1
Å0
*
( p / 2 ) Du p / u
p (cos u1 0 cos u2 ) 2 p iDu 2
*
u2
u1
sin u du. u
[3]
Now p takes i 0 1 values from 1 to i 0 1 as there are only i 0 1 down triangles in the band. Evidently, d v, is also independent of the index p. The total solid angle subtended by the given band is equal to id vD / (i 0 1)d vÇ Å (p /2)(cos u1 0 cos u2 ), as expected. Both d vD and d v, involve the evaluation of the integral (sin u / u ) du. This integral can be expanded into a fast converging series so that it can be easily evaluated numerically. However, for typical partition numbers (N) used in the simulation of randomly oriented powder spectra, the difference between d vD and d v, is very small. For example, when N Å 19 ( Du Å 57 ) dv, 0 d vD/d v, is 1.76 1 10 02 for the band defined by u1 Å 857 and u2 Å 907, the band in which the largest difference is found. This difference diminishes to 1.27 1 10 03 for the band defined by u1 Å 57 and u2 Å 107, the band of the smallest difference. With a larger N, the difference becomes smaller, which is as expected. Taking N Å 91 as another example, ( Du Å 17 ), the respective values are 3.7 1 10 03 and 5.08 1 10 05 . Each triangle in Fig. 2a can be easily subpartitioned into smaller triangles, referred to as tiny triangles. In Fig. 2b, a selected triangle is further partitioned into 81 tiny triangles with a subpartion number M Å 10. The grid formed in such a subpartition can still be described by Eq. [1]. In this particular case, C1 is stepped in a smaller step of p /[2(N 0 1)(M 0 1)] from C1 Å 457 to C1 Å 547, the two corresponding curves which bound the triangle. A similar process is applied to curves in sets 2 and 3. The subpartition is used in conjunction with a linear interpolation scheme. Although we can calculate the solid angles corresponding to each of these tiny triangles by using Eq. [2] and Eq. [3], such exact calculations are hardly justifiable in view of the approximations involved in the linear interpolation. For example, for the N Å 19 partition with a subpartition number M Å 51 (in other words, each triangle is subpartitioned into 2500 tiny triangles), the biggest difference among the solid angles subtended by the tiny triangles in a given triangle was found to be less than 5%. Therefore, in our program, we assumed all the tiny triangles in a given triangle subtend the same solid angle, thus saving many calculations which would be required if Eq. [2] and Eq. [3] were used.
AP: Mag Res
4
WANG AND HANSON
A NEW INTERPOLATION SCHEME
Most interpolation schemes adopted so far are two-dimensional (6, 7, 15) with u and f as variables. One-dimensional interpolation has been considered by Mombourquette and Weil (10) but they only considered the simplest interpolation method, linear interpolation. In most cases, linear interpolation cannot satisfactorily reconstruct the experimental spectra without performing exact calculations of the positions and intensity at a very large number of orientations. When higher-order interpolation methods are considered, first derivatives with respect to u and f are often required. The task of calculating these first derivatives can easily become as formidable as obtaining the position/intensity information. An important feature in our new interpolation scheme is that we only use one-dimensional interpolation methods in which the demand for the first-derivative calculations is dramatically reduced. Specifically, we use cubic spline interpolation for data points along the edges of the triangles and linear interpolation for those inside the triangles from the data obtained from the cubic spline interpolations. Global Cubic Spline Interpolation Cubic spline interpolation is featured by its smoothness in the first derivative and continuity in the second derivative, both within an interval and at its boundaries. Therefore, it should be a good interpolation method for the data points along the edges of the triangles. Strictly speaking, among the three sets of curves in SOPHE, only set 1 (Eq. [1a]) involves one variable, f. The other two sets of curves involve both u and f. However, by using Eq. [1b] and Eq. [1c], sets 2 and 3 can also be treated as involving only one variable. In both sets, we keep u as the variable and treat f as a parameter and thus one-dimensional interpolation methods can be used for all three sets. The cubic spline interpolation process is divided into two steps (16). For a given set of data, xi and yi Å y(xi ), i Å 1, . . ., n, and two end conditions, the first step is to calculate the second derivatives d 2 y/dx 2 (or the first derivatives, dy/ dx) at the points xi (i Å 1, 2, . . ., n). Different types of end conditions can be used (16) and here we adopt the firstderivative end condition, which is to specify the two first derivatives at the two endpoints x1 and xn . The second derivatives are then obtained with a C function called spline in ‘‘Numerical Recipes in C’’ (17). For a given value x, the second step is to find the interval within which the specified x is located. In the algorithm we developed, this location search has been avoided, thus saving some computing time. In the actual calculation, we store the second derivatives. There are 6(N 0 1) first derivatives, 2(N 0 1) for each set of interpolation. This number is significantly reduced from what is required in the two-dimensional schemes used by
m3839$0679
10-10-95 22:18:52
magas
previous authors (6, 7). For example, in the interpolation scheme used by Gribnau, van Tits, and Reijerse (7), a total of 3N 2 first derivatives were required. For a partition number of N Å 31, it requires 16 times as many first-derivative calculations as used in our scheme. For each set, there are N 0 1 cubic spline interpolations with the number of data points from 2 to N in steps of 1. We have tested the quality of the interpolated data from the cubic spline interpolation by comparing them with those obtained through exact calculations on a number of spin systems and the results were found to be excellent. Here, we examine a specific example of an S Å 32 spin system. The spin Hamiltonian used has the following form K Å g bBrS / D [S 2z 0 13 S(S / 1)] / E(S 2x 0 S 2y ),
[4]
where the terms have their usual meanings (2). The comparison was made on a highly anisotropic system in which E/D Å 13, the so-called fully rhombic symmetry. We used the data for a Cr(CN) 30 center in a KCl lattice (18). For 6 this center, D Å 2.73 GHz, g à 1.990, and E/D Å 13. The microwave frequency was set to 9.5 GHz. We made the comparison on the three allowed transitions, 1 } 2, 2 } 3, 3 } 4 with the energy levels numbered in order of increasing energy. In the exact calculation, we used a bisection method in locating the resonant positions with an accuracy of {0.01 mT. It should be pointed out that the bisection method is one of the least efficient methods in locating resonant fields and is used simply for the purpose of demonstrating the SOPHE algorithm. The resonant field position of the allowed transitions varies between 160 and 536 mT. Table 1 presents a summary of the comparison. It also includes the maximum and minimum resonant fields for each transition. The data points selected for comparison were obtained by stepping u by 0.57 for curve sets 2 and 3, while points of separation with approximately equal arc lengths of 0.57 were selected for curve set 1. The overall agreement in Table 1 is excellent. Naturally, when N gets larger, the grid becomes finer and the agreement is better. With N Å 31, no data points with differences between the calculated and interpolated larger than 0.02 mT were found for transition 3 } 4 and only 12 such points were detected for transitions 1 } 2 and 2 } 3, which represents only 0.16% of the total number of points compared. The average deviation is very small, a few microtesla for all three transitions. A similar comparison can be made for the intensity but normally the transition probability varies much less than the resonant field positions and so this was not included in Table 1. Another important feature worth mentioning explicitly about Table 1 is that normally one would expect a better agreement for transition 2 } 3 as its resonant field varies
AP: Mag Res
5
SIMULATING RANDOMLY ORIENTED POWDER SPECTRA
TABLE 1 Comparison of the Resonant Fields Interpolated through the Global Cubic Spline Interpolation (Binterpol) with Those Calculated by Full Matrix Diagonalization (Bexact) N
a
dub (degrees)
Npc
3}4
Bmaxg Å 536.0039 (y axis)
19 25 31
0.01 0.01 0.01
1}2
Bmax Å 536.0039 (z axis)
19 25 31
0.01 0.01 0.01
2}3
Bmax Å 330.5821 (y or z axis)
19 25 31
0.01 0.01 0.01
4842 6488 7480
4842 6488 7480
4842 6488 7480
K1d
K2e
Average deviation f (mT)
Local Linear Interpolation
Bmin Å 160.7016 (z axis) 0 0 0
3 0 0
0.90 0.31 0.16
Bmin Å 160.7016 (y axis) 19 9 0
71 26 12
1.97 0.69 0.33
Bmin Å 299.7798 (x axis) 19 7 0
70 25 12
1.68 0.59 0.29
a
N is the partition number. du is the step used for evaluating the first derivatives. c Np is the total number of data points compared. d K1 is the number of data points for which DB § 50 mT (DB Å ÉBinterpol 0 BexactÉ). e K2 is the number of data points for which DB § 20 mT. f Average deviation is defined as the average of DB for all the points compared. g Bmax (mT) and Bmin (mT) are the maximum and minimum resonant field positions when the magnetic field is applied along the axes shown in brackets. b
only in a small fraction of the field range spanned by the other two transitions, but the opposite is shown when comparing transitions 2 } 3 and 1 } 2. The reason for this is that the quality of the interpolated data depends on the accuracy of the first derivatives. As said, the first derivatives in the full diagonalization method must be obtained numerically, and to evaluate all the first derivatives to a good precision is not an easy task. This is because first derivatives normally take values over a wide range. In our test case, they vary from a few hundred microtesla/degree to over 10 mT/degree, and as they are not known a priori, it is impossible to choose a good value for du or df for each firstderivative evaluation. In our program, we set the same value for du for all first-derivative calculations. For curve sets 2 and 3, an appropriate amount of df is necessary, which is obtained through Eq. [1b] and Eq. [1c]. For curve set 1, we choose values for df by using df Å du/sin u. We have
m3839$0679
10-10-95 22:18:52
tested the quality for the interpolated data as a function of du for all three transitions and the results for transition 1 } 2 are shown in Table 2. In general, any values between 0.0057 and 0.17 are acceptable. However, we found 0.017 is an optimal value for du for our test case with the resonant field positions determined to an accuracy of {0.01 mT.
magas
After the global interpolation, the position/intensity/linewidth information can be viewed as a ‘‘skeleton’’ based on the designed grid in Fig. 2a. Further subpartitioning can be made for each triangle, and for any data point inside the triangle, we employ a linear interpolation method. This process can be viewed as a ‘‘tile-filling’’ process. One advantage of the linear interpolation method is that it does not involve any derivatives, so that no further exact calculations are required. The justification of the linear interpolation relies on the anisotropy of the system under consideration and such a justification should always be achievable by varying the grid size. A number of linear interpolation schemes can be employed, which are schematically shown in Fig. 3. The singleset scheme shown in Fig. 3a suits most situations where the resonant field position hardly changes along one of the edges. In this case, we linearly interpolate the data points with the set which includes the edge along which the resonant field position shows little change. For other situations, the combined (Fig. 3b) or the averaged scheme (Fig. 3c) can be used. Obviously, the averaged scheme requires more computing time. However, we noticed that the averaged scheme in general produced higher-quality (less ripple) simulations
TABLE 2 Comparison of the Resonant Fields Interpolated through the Global Cubic Spline Interpolation Method (Binterpol) with Those Calculated by the Full-Matrix-Diagonalization Method (Bexact) for Transition 1}2 with Different Values of dua N
du (7)
Np
K1
K2
Average deviation (mT)
31 31 31 31 31 31 31
0.005 0.01 0.05 0.1 0.5 1.0 3.0
41,290 41,290 41,290 41,290 41,290 41,290 41,290
0 0 0 7 496 2632 6902
57 55 66 116 2091 4873 10058
0.33 0.31 0.51 0.82 3.50 53.3 98.0
a The symbols have the same definitions as in Table 1. Here, a smaller step of 0.17 (0.57 was used for Table 1) was used in selecting the points for comparison in order to show there was absolutely no oscillations in the interpolation. Consequently, the total number of data points increases from 7480 to 41,290 for N Å 31. The vertex points were excluded in all comparisons.
AP: Mag Res
6
WANG AND HANSON
FIG. 3. Schematic representations of the possible linear interpolation schemes. (a) Single-set scheme. (b) Combined scheme. (c) Averaged scheme.
than the combined scheme. Thus, in our program, the averaged scheme is used. In the averaged scheme, each point inside the triangle is interpolated three times from three different sets and an averaged value is taken for each point. RECONSTRUCTION OF A RANDOMLY ORIENTED POWDER SPECTRUM
In this section, we present some details of the computer implementation of SOPHE. Suppose the position/intensity/ linewidth information for a selected transition has been obtained at the vertex points for a given partition number N through a certain method. We can then start to calculate the first derivatives required for the cubic spline interpolation. If only the resonant position is globally interpolated, there are 6(N 0 1) first derivatives, 2(N 0 1) for each set of global interpolations. If both resonant position and intensity are globally interpolated, the number of first derivatives is simply doubled. Second derivatives are obtained with the function spline (17). These second derivatives, together with the position/intensity/linewidth information at the vertex points, are passed on to a function which handles interpolation and convolution. The interpolation process is made triangle by triangle and band through band. It starts with band i Å 1 at the top of the grid and progresses toward the xy plane. For a given band, the data elements corresponding to the vertex points can be arranged as a one-dimensional array
m3839$0679
10-10-95 22:18:52
magas
so that to move from one triangle to the next is simply to shift the index by 1. Figure 4 shows such an indexing scheme for a band with five triangles. The three vertex points which define a triangle are j, j / 1, j / 2 with j shifting from 1 to 5. Therefore, the whole interpolation/convolution process involves only repeating the same routine which handles a single triangle and that is how we implemented the SOPHE method. In the actual calculation, we first obtain a profile of wkl d v against the resonant position for a selected triangle; wkl represents the intensity for the transition labeled by kl. This is done by going through all the tiny triangles into which the given triangle is partitioned. Then the profile obtained is convoluted with a given lineshape with the linewidth averaged from those at the vertex points before adding its contribution to the final spectrum. This process is repeated for every transition. Since EPR is a field-swept technique, corrections for the anisotropic transition probability and linewidth must be dealt with in a proper way in computer simulations. This is especially important for high-spin systems. When perturbation theory is used, these corrections can be taken into account point by point. This normally involves calculation of the spectrum in frequency space and then transforming back into field space and requires very little additional computational time (2). However, in the full-matrix-diagonalization method, such point-by-point corrections become impractical. We have, instead, used an approach introduced by van Veen (6) in which the corrections for the anisotropic transition probability and linewidth are applied at the triangle level, i.e., the vertex points where exact calculations are performed. An S Å 52 spin system with g Å 1.990, D Å 1.0 GHz and E/D Å 13 is used to demonstrate the efficiency of this algorithm based on the SOPHE method. Here, again we assume full rhombic symmetry with a large zero-field splitting, 7.06 GHz, which is comparable with the microwave frequency used, 9.0 GHz. A partition number N Å 19 was used, which means exact calculations were performed at 190 vertex points. In Fig. 5a, the interpolation was made in steps of 17, corresponding to 3905 interpolated points per transition for each of the five allowed transitions. As shown, the reconstructed spectrum is very ‘‘noisy.’’ In Fig. 5b, the interpolation was made in steps of 0.17, resulting in 405,260 data
FIG. 4. The indexing scheme used in our computer implementation of the SOPHE method.
AP: Mag Res
SIMULATING RANDOMLY ORIENTED POWDER SPECTRA
7
points per transition and the computer ‘‘noise’’ has been removed nearly completely. Figure 5b consumed 120 s more, 24 s/transition, in CPU time on our Sun SPARCstation 10/ 30 workstation. This means that at a marginal cost of 24 s, more than 400,000 data points have been generated. In fact, the large proportion of the consumed CPU time was spent in the sorting routine and the data-generating speed in the interpolation routine was about 50,000 data points per second, a remarkable achievement for the SOPHE method. Further improvement in spectrum quality can be made by using a finer grid (a larger value for N). In Fig. 5c, a partition number of N Å 31 ( Du Å 37 ) was used with a subpartition of 0.17. DISCUSSION AND CONCLUSION
5
FIG. 5. Simulated spectra for a model S Å 2 spin system (g Å 1.990, 1 D Å 1.0 GHz, and E/D Å 3). (a) Parameters used in the simulation were partition number, N Å 19 (190 vertex points); subpartition number, 6 (3905 interpolated points per transition); field axis resolution, 4096 points; field sweep-width, 350 mT; center field, 325 mT; microwave frequency, 9.0 GHz. A Gaussian lineshape with a half-width at half-maximum of 30 MHz was used. The linewidth was assumed to be isotropic. (b) Parameters used in the simulation were the same as for (a) except for a larger subpartition number of 51 (405,260 interpolated points per transition). (c) Parameters used in the simulation were the same as for (b) except with a finer grid, N Å 31 (496 vertex points).
m3839$0679
10-10-95 22:18:52
magas
The most significant features of this new algorithm based on the SOPHE method are its efficiency in generating a very large number of data points from those at vertexes, which are normally obtained through a time-consuming process, and its reliability. Data can be interpolated extremely rapidly, tens of thousands of data points per second on modern computers. For example, on our Sun SPARCstation 10/30, the data-generating speed in the interpolation process is Ç50,000 points/s. Thus, the SOPHE method has reduced the computing task mainly to calculating the position/intensity information at the vertex points, computing the first derivatives, sorting, and convoluting. This is a remarkable improvement in computer reconstruction of randomly oriented powder spectra. In particular, the SOPHE method should significantly ease the burden of simulating randomly oriented powder spectra of spin systems of large dimensions. The quality of the interpolated data has been demonstrated to be very high. The interpolated values on the edges can be expected to be very close to the true values if the first derivatives are computed with a good precision as typified in Table 1. This is important as the globally interpolated data are to form the base on which local linear interpolations are performed. Although we have only discussed the application of the SOPHE method to orthorhombic spin systems, it can be easily extended to monoclinic (C2 , C2h , and Cs ) and triclinic (C1 ) symmetries where two and four octants are required respectively. In fact, the SOPHE method introduced in the previous section is a special case of the more general situation where a section of the unit sphere is defined as u1 Å 07 and u2 Å 907, and f1 and f2 are two arbitrary angles with f2 ú f1 . Equation [1] can still be used providing f is replaced by f 0 f1 and p /2 is replaced by f2 0 f1 in both Eq. [1c] and the constant C2 . Therefore, for monoclinic symmetries, the SOPHE method can be applied by choosing f1 Å 07 and f2 Å 1807 if the z axis is defined as the symmetry axis (C2 and C2h ) or as the axis perpendicular to the reflection
AP: Mag Res
8
WANG AND HANSON
plane (Cs ). In the case of triclinic symmetry, half a sphere must be considered, which means f1 Å 07 and f2 Å 3607. In this case, set 1 curves are closed curves and periodic end conditions must be used. Alternatively, in monoclinic and triclinic symmetries, each individual octant can be treated separately in the same manner as for orthorhombic symmetry and their contributions are then added to the final spectrum. Subpartitioning of each triangle can also be easily made on an individual basis; that is, how fine each individual triangle is partitioned into tiny triangles is determined by the difference in the resonant field positions at the three vertex points. In other words, it can be determined by the ‘‘local’’ anisotropy for a given transition. Thus, the SOPHE method has a great flexibility in a sense that it handles the anisotropy step by step through the space, and, of course, transition by transition. The SOPHE method is a general method for simulating randomly oriented magnetic resonance spectra. It can be used for any computer reconstruction where the spectral function must be sampled at a very fine grid of the unit sphere, in other words, at a very large number of data points, and sampling at each point takes a considerable amount of computing time. For example, the method should find its use in simulating randomly oriented powder spectra in solid state NMR, electron nuclear double resonance (ENDOR), electron spin-echo envelope modulation (ESEEM), and nuclear quadrupole resonance in solids. Specific details concerning the availability of the SOPHE software can be obtained by writing to the authors.
REFERENCES 1. P. C. Taylor, J. F. Baugher, and H. M. Kriz, Chem. Rev. 75, 203 (1974), and references therein. 2. J. R. Pilbrow, ‘‘Transition Ion Electron Paramagnetic Resonance,’’ Clarendon Press, Oxford, 1990, and references therein. 3. G. M. Muha, J. Magn. Reson. 49, 431 (1982); 53, 85 (1983). 4. Y. Siderer and Z. Luz, J. Magn. Reson. 37, 449 (1980). 5. L. Andreozzi, M. Giordano, and D. Leporini, J. Magn. Reson. A 104, 166 (1993). 6. G. van Veen, J. Magn. Reson. 38, 91 (1978). 7. M. C. M. Gribnau, J. L. C. van Tits, and E. J. Reijerse, J. Magn. Reson. 90, 474 (1990). 8. G. G. Belford, R. L. Belford, and J. F. Burkhalter, J. Magn. Reson. 11, 251 (1973). 9. A. Kretter and J. Huttermann, J. Magn. Reson. 93, 12 (1991). 10. M. J. Mombourquette and J. A. Weil, J. Magn. Reson. 99, 37 (1992). 11. L. Cugunov, A. Mednis, and J. Kliava, J. Magn. Reson. A 106, 153 (1994). 12. M. Krusche, M. Neu, and H. Ziegler, J. Magn. Reson. 95, 368 (1991). 13. E. Anderson, Z. Bai, C. Bischof, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, S. Ostrouchov, and D. Sorensen, ‘‘LAPACK Users’ Guide,’’ Society for Industrial and Applied Mathematics, Philadelphia, 1992. 14. (a) M. J. Nilges, Ph.D. thesis, University of Illinois, Urbana, Illinois, 1979; (b) R. L. Belford and M. J. Nilges, Computer Simulation of Powder Spectra, EPR Symposium 21st Rocky Mountain Conference, Denver, Colorado, 1979; (c) A. M. Maurice, Ph.D. thesis, University of Illinois, Urbana, Illinois. 15. D. W. Alderman, M. S. Solum, and D. M. Grant, J. Chem. Phys. 84, 3717 (1986).
ACKNOWLEDGMENTS
16. Su Bu-qing and Liu Ding-yuan, ‘‘Computational Geometry— Curves and Surface Modeling,’’ Academic Press, San Diego, 1989.
We thank Doug Maclachlan and Kevin Gates for providing us some fast diagonalization routines and useful references, and Simon Teed for drawing some of the figures. The Australian Research Council is also thanked for awarding Research Grant A29331095 to G.R.H. which made this work possible.
17. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, ‘‘Numerical Recipes in C,’’ p. 113, Cambridge Univ. Press, Cambridge, 1992.
m3839$0679
10-10-95 22:18:52
magas
18. D. M. Wang, D. H. Hutton, and J. R. Pilbrow, J. Phys. C, 19, 789 (1986).
AP: Mag Res