Determination of the g-, hyperﬁne coupling- and zero-ﬁeld splitting tensors in EPR and ENDOR using extended Matlab codes

The analysis of single crystal electron magnetic resonance (EMR) data has traditionally been performed using software in programming languages that are difﬁcult to update, are not easily available, or are obsolete. By using a modern script-language with tools for the analysis and graphical display of the data, three MatLab (cid:1) codes were prepared to compute the g, zero-ﬁeld splitting ( zfs ) and hyperﬁne coupling ( hfc ) tensors from roadmaps obtained by EPR or ENDOR measurements in three crystal planes. Schonland’s original method was used to compute the g- and hfc -tensors by a least-squares ﬁt to the experimental data in each plane. The modiﬁcations required for the analysis of the zfs of radical pairs with S = 1 were accounted for. A non-linear ﬁt was employed in a second code to obtain the hfc -tensor from EPR measurements, taking the nuclear Zeeman interaction of an I = ½ nucleus into account. A previously developed method to calculate the g- and hfc -tensors by a simultaneous linear ﬁt to all data was used in the third code. The validity of the methods was examined by comparison with results obtained experi- mentally, and by roadmaps computed by exact diagonalization. The probable errors were estimated using functions for regression analysis available in MatLab. The software will be published at https://doi.org/10. 17632/ps24sw95gz.1, Input and output examples presented in this work can also be downloaded from https://old.liu.se/simarc/downloads?l=en. (cid:3) 2021 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
The structure and reactions of paramagnetic defects in solids have been extensively examined by electron magnetic resonance (EMR) methods for more than 60 years. Methods to extract the g-tensor by regular electron paramagnetic resonance (EPR) measurements on single crystals were presented in early work [1], and further developed in software containing error analyses and modifications to treat electron-nuclear double resonance (ENDOR) data [2][3][4][5]. However, these codes may no longer be available or maintained. Sophisticated software packages have since then been developed to aid the interpretation of EMR data [6,7], but procedures to extract the coupling tensors from single crystal experiments do not seem to be available. The short MatLab Ò codes presented here do for simplicity or practical reasons neglect some options accounted for in more extensive software [2][3][4][5], like to include data obtained at different spectrometer frequencies, or at different electronic states, e. g. at m S =±½ in ENDOR measurements of free radicals. Data from the low frequency branch can be analysed if required but the intensity was often low, and only the data of the high frequency branch could be employed in the present work. To improve the accuracy of the g-and hfc-tensors a slightly simplified method was implemented for the correction of errors due to crystal misorientation than those considered in [3,5]. The programs presented here may be of interest in applications concerned with radiation research, biochemical and biophysics applications and other fields where mainly free radicals and triplet state molecules in organic substances are observed. In the following and in the software we therefore limit the treatment to cases with electron spin S = ½ and 1. The spin-Hamiltonian Eq. (1) with electron and nuclear Zeeman terms, an hfc term and for S = 1 a zfs term S. D.S, containing only quadratic terms, has then been employed [8 pp 225-253], see also [9 pp 197-199].
H ¼ l B B:g:S þ S:D:S þ S:A:I À l N g N B Á I The determination of the g and A (in fact g 2 and A 2 ) and D tensors in general requires EPR and/or ENDOR measurements with the magnetic field vector B in a set of 2-4 planes followed by fitting of the data to ellipsoidal surfaces. The selection of these planes depends on the specific samples available, their crystal symmetry https://doi.org/10.1016/j.jmr.2021.106956 1090-7807/Ó 2021 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). and the local site symmetry (LSS) determined by the specific defect or radical geometry in the crystal. Obviously the LSS is usually lower than the crystal symmetry. Often these factors will naturally lead an experienced experimentalist to the choice of a laboratory axis system (LAS) labeled xyz to describe the measurement data (EPR magnetic fields and/or ENDOR frequencies). Examples will be discussed below. The raw experimental data depend on one or more of the symmetric second rank tensors mentioned above [10,11]. When the second rank properties can be determined for a sufficient number of orientations of the magnetic field, the transition from the LAS to the principal axis system (PAS) of the tensor involved is mathematically straightforward (reduction of an ellipsoid) and will lead to diagonalization and knowledge of the principal values and the corresponding principal axes. In the simple case where S = ½ and where hfc interactions are absent (I = 0) or negligible, the g-values can be simply calculated from the EPR experimental values and the complete g-tensor can be easily determined, of course, with some experimental uncertainty. In the next section, a possible practical procedure is described in more detail (Schonland).
When hfc and/or zfs interactions become important, the analysis of the experimental data, influenced by more than just the gtensor, becomes more difficult [11,12]. Historically a number of practical methods have been developed and the spectroscopist has to decide which one to use, taking into account the relative strengths of the interactions, the theoretical approximations used, also in comparison with the experimental data available and their unavoidable experimental errors.
In the LAS the S.D.S term in eq. (1) can be expanded in a series of diagonal and off-diagonal terms or, alternatively, the H ZFS can be written as [11,16,17]: In the PAS of the D-tensor, the following equivalent forms have been reported in terms of the Extended Stevens operators : The b q k are scaled parameters that are also often used [11,13]. In cases of tetragonal, hexagonal or trigonal LSS symmetry, the D-tensor becomes axial and only one zfs-parameter is required (D or B 0 2 ). When the LSS symmetry is orthorhombic or lower, it is advisable to take the standardization into account as recommended by Rudowicz and Bramley [14], see also [13]. This allows a correct comparison between data. To keep the ratio k = E/D in the range 0-1/3, the D-tensor principal values should be ordered as follows: jD z j ! jD x j ! jD y j Important in this respect is the fact that, if more than one tensor has to be taken into account (e.g. g and A), the PASs do not necessarily coincide anymore. E.g. in cases of monoclinic or triclinic local symmetry, low symmetry effects may occur [15].
In the original (and present) software and for the radical systems mainly envisaged, situations with non-coincident tensors were not included. The software was less intended for high spin systems and therefore S was considered to be not larger than 1. The influence of the g-anisotropy on the values of the hfc -and zfs -tensors discussed in Section 5 was systematically assumed to be negligible, avoiding problems with non-coincident PAS and more complex Hamiltonians.
Moreover, it has to be considered that further theoretical approximations were made, including 1st or 2nd order perturbation formulas, allowing to access g, A and D-values (via experimental splittings) more easily (but less accurately). A nice demonstration of (the advantage of) perturbation theory is the determination of the g-values from central line positions when hfc interaction is present (see SAR-example below). Thereafter the determination of the g-tensor is again straightforward.
Also for the A and D-tensors, this leads to mathematical treatments that are very similar or even identical to single (g) tensor situations. Examples will also be given in the next section.
Historically introduced by lack of fast computers, the approximate formulas provide a lot of insight in the experimental data patterns to be analyzed, and it has to be borne in mind that the errors introduced are often smaller than the experimental error. E.g., the 2nd order correction for a hfc coupling of about 20 MHz was estimated to be of the order 0.01 MHz [8 pp 193-194]. An analysis of the ENDOR data of a weakly interacting proton and of the relatively large H a splitting typical of an aliphatic radical has indicated for example ( Table 5 in this work) that the errors caused by the use of 1st order perturbation theory were of the order 0.05 and 0.2 MHz, respectively. The smaller values are considered as acceptable errors, while the larger deviation might be reduced by modelling the resonance frequencies by 2nd order perturbation theory [23] as suggested in Section 5. Of course, the potential users of the software have to decide whether the advantages (ease of use, direct availability, transparency, insight, . . .) cancel or even outweigh the disadvantages.
A further development of the software based on the more general Hamiltonians used e. g. in [6,7] is desirable to establish a common standard for the determination of tensor quantities in single crystal EMR studies. This subject is not within the scope of the present work, however.
In the next section, a first possible practical procedure is described in more detail.

The Schonland method
The Schonland method to obtain the g-tensor [1] involved EPR measurements in three crystal planes, by rotating the crystal about an axis perpendicular to the magnetic field. The crystal symmetries that are considered are orthorhombic, monoclinic or trigonal but the LSS will in most cases lead to a triclinic g 2 and A 2 tensor.
A set of measurements obtained by rotation of the magnetic field in a certain plane can be represented by the intersection of this plane with the tensor representation ellipsoid, i.e. an ellipse. The latter can be easily fitted with three parameters, the maximum and minimum g-values in the plane and the angle at which the maximum g-value occurs [1]. When the rotation plane contains a principal axis, the corresponding principal value can be directly determined. In general three sets of measurements (three ellipses) will allow the reconstruction of the ellipsoid and thus the g-tensor. When g anisotropy is neglected, as will be assumed in the present paper, also the A and D-tensors can be determined in the same way, using formulas requiring no further approximations.
The (x,y,z) LAS that we used for the measurement of orthorhombic crystals is the natural orthogonal (a, b, c) system (crystal axes system, denoted by CAS). For the crystals with a monoclinic space group, the (a, b, c*) or (a*, b, c) axes, with c* or a* perpendicular to b and to a or c (denoted by CAS*) were chosen to analyze the experimental data, as will be demonstrated below.
As explained above, the g-factor depends on the direction of the magnetic field according to the formula in Eq. (5) or an equivalent one in terms of the maximum and minimum g-values in the plane and the angle at which the maximum g-value occurs [1]. The orientation of the magnetic field B in each plane is specified by the angle h.
The maximum and minimum values g 2 AE of this expression are given by Eq. (6) and are found for values of h that are 90°apart.
For an orthorhombic crystal, the elements of the tensor T = g 2 were expressed in the CAS (x, y, z) = (a, b, c) and calculated by a least-squares fit of Eq. (5) to the data in the bc, ac and ab crystal planes respectively. The principal g-values were calculated as the square root of the corresponding g 2 values. As is well-known, the principal directions of g 2 and g coincide [1]. For a monoclinic crystal, measurements in the ab, ac, and bc planes were employed with the tensor elements specified in the CAS* (a, b, c*) system, Fig. 1 [1]. The procedure to evaluate the tensor elements from measurements in crystal planes (bc*, ac and a*b) perpendicular to the a, b and c crystal axes, i.e. the three non-orthogonal rotation axes with orientations as in Fig. 1, is analogous. The tensor elements were then eventually expressed in the orthogonal LAS (x, y, z) = (a*, b, c). All tensor elements except T yz were obtained from rotations about the b and c axes. The quantity T 0 yz was measured in the bc* plane and expressed in terms of the already determined elements, Eq. (7), taking in account the transformation properties of Cartesian tensors. The T zz value in Eq. (8) was obtained in an analogous manner and employed as a consistency check. The angle e = b À 90°was introduced in Eqs. (7) and (8) to emphasize the similarity with the equations in [1].
The orthorhombic case corresponds with the case e = 0.
The procedure was verified using simulated calculated g-tensor roadmaps. The fits obtained with constant and adjustable weights of the data with the Matlab functions regress and robustfit then overlapped. The principal values and direction cosines used to generate the roadmap data in Fig. 1 were recovered with negligible error from the rotations about the a, b and c axes of a monoclinic crystal with an angle b = 110°between the a and c axes.
Performance: The data in [18] were reanalyzed in order to test the performance on the experimental g-factor data of the stable alanine radical (SAR) marked with o in Fig. 2.
The g-tensor in Table 1a, obtained using regress fits in rotations about the a, b and c orthorhombic axes, deviated somewhat from the reported ones. Good agreement with Ref. [18] was obtained by reversing the sense of rotation about c, Table 1b. In this and tables below the tensor orientation with respect to the LAS was specified with the eigenvectors for each principal value. The corresponding Euler angles, displayed on the computer screen at runtime, can also be stored if required. . The difference is probably due to the use of first order theory to obtain the ENDOR frequencies in the Schonland function. The tensor obtained by analyzing the roadmaps of the weaker coupling generated by 1st order theory agreed to within 0.01 MHz for the principal values and to at least 3 decimal places for the direction cosines. The higher order corrections for the weaker coupling were assumed to be insignificant. The procedure to extract the hfc tensors of organic radicals with the Schonland method by analysis of the ENDOR frequencies was extended in [5,8 Ch. 4] and works cited therein. Simple expressions for the hfc tensor elements were obtained by employing data for the high as well as the low frequency branch. Errors due to the uncertainty of the frequencies and the crystal orientation were both accounted for. The Schonland function used in this work was designed to obtain the hfc tensor by analysis of the high frequency branch, given that the signal intensity of the low frequency branch was usually weak. Data from the low frequency branch can, however, be analyzed if required, as verified using simulated roadmaps for the analysis. The probable errors were estimated using functions for regression analysis available in MatLab.
The performance for the analysis of ENDOR data was investigated by reanalyzing the data of an X-irradiated ammonium tartrate crystal. The measurements were performed by rotation of the crystal about the a, b and c crystal axes, with monoclinic angle b = 92.4° [20]. Fits with constant and adjustable weights of the data were obtained using the Matlab functions regress and robustfit. The fits were nearly identical, Fig. 3. The principal values and directions of the R1(3) tensor obtained by the Schonland function in Table 2 coincided exactly with those computed with the Nelson function in Table 8, indicating that the two methods were consistent provided that the orientation of the magnetic field with respect to the crystal axes is accurately known. The automatic method to correct for crystal misalignment described in Sections 2 and 3 is not applicable with the Schonland method, however, and the angle adjustments calculated by the Nelson function (Table 8) were therefore inserted manually. Note also that the principal values of the hfctensor in Table 2 were obtained with reversed signs to those reported in [20] for the R1(3) tensor. This is not unreasonable considering that the absolute signs were not determined by the analysis. The Euler angles specifying the orientation of the fitted tensor in Table 2, (-5.76 50.16 -131.34)°, were nevertheless in good agreement with the experimental values, (-5.53 48.39 -131.41)° [ 20] (see the actual data for the R1(3) tensor presented in Table 5).
The performance of the Schonland function for the analysis of the zero-field-splitting (zfs) of a radical pair was investigated by employing the data of an X-irradiated potassium dithionate crystal [21]. The X-and Q-band spectra observed by rotation about an axis (x) perpendicular to the trigonal axis of the crystal contained a weak signal assigned to a (SO 3 À ) 2 radical pair with an axially symmetric zfs, Fig. 4. The analysis was made with a slight modification of the Schonland method, by fitting of the a, b and c values (see formula (2)) to the splitting between the m S = 1?0 and 0?À1 transitions rather than to its square. Performance test of the Schonland function for the analysis of the g-factor variation in three orthogonal planes (o) of the SAR radical in an irradiated L-alanine crystal [18]. The fits in solid lines obtained with constant and adjustable weights of the data employing the Matlab functions regress and robustfit did nearly coincide. Fits were obtained using the Matlab functions regress with equal weights of the data, and robustfit, which gives lower weight to points that do not fit well by an iteratively reweighted least squares algorithm. Assuming that the axis of the radical pair coincides with the trigonal crystal axis, axial D-tensor symmetry was assumed [21]. In axial symmetry the smallest splitting always cor-responds to -D. The perpendicular component can in principle be directly obtained from the data in Fig. 4, but due to the restricted data set and the substantial extrapolation, only an approximate value of D=|F -|= 3.23 mT was estimated using an expression analogous to Eq. (6) to calculate the maximum and minimum zfs, F ± .
The hfc EPR data of organic free radicals can in many cases be analyzed by neglecting the influence of the nuclear Zeeman interaction. There are, however, exceptions where the term is not negligible, like for the proton couplings in p-electron radicals; another case is radicals containing fluorine atoms. An outer and inner doublet of EPR lines due to the hfc of 1 H or 19 F nuclei (I = ½) were observed in these cases. A non-linear fitting procedure (EPR6) to the outer doublet splitting was developed to obtain the a-H proton tensor of the free radical formed by H-atom abstraction in cirradiated single crystals of barbituric acid monohydrate [2]. We have previously employed the theory derived by Weil and Anderson [22] to extract the principal values and orientations of the hyperfine coupling tensors of nuclei with I = ½ in an analysis of the EPR spectra of irradiated crystals [4]. A similar treatment was adopted here, using the formulae due to Iwazaki [23, 8 pp 175-181] in an arbitrary coordinate system to compute the separation,  [20] in rotations about the a, b and c monoclinic axes. The fits in solid lines, obtained with constant and adjustable weights of the data employing the Matlab functions regress and robustfit, overlapped closely.  4. Roadmap of the zfs of an (SO 3 À ) 2 radical pair in an irradiated potassium dithionate crystal [21]. The solid lines were fitted using the MatLab functions regress and robustfit, to the zfs-data marked with circles,1 G = 0.1 mT.

S (mT), of the outer doublet lines.
A is the hfc -tensor, b = -g N Ál N ÁB/ (gÁl B ) the nuclear Zeeman splitting, both in mT units, V a unit tensor and l a unit vector oriented along the magnetic field.
Eq. (9) was obtained by neglecting the g-anisotropy in [23]. The direction cosines employed in the fits were obtained from Eq. (10), by adapting a formula in classical mechanics to compute the new position of a point after rotation about a fixed axis N [24].
The direction of l, upon rotation over an angle h about N, was measured relative to an initial reference direction l 0 of the magnetic field. Following the convention in [1], I 0 is a unit vector directed along x, y and z, respectively, upon rotation about the z, x, and y axes. The angle h 0 corresponds to l||l 0 . The corresponding reference direction for an arbitrary crystal plane must be provided by its direction cosines.
The method was tested by analysis of the hfc splittings due to the a-H proton of radical III in irradiated L-asparagine. Simulated EPR roadmaps were calculated at a microwave frequency of 9.7 GHz, using the hfc tensor reported in [19], either by adopting the first order treatment in [22], or the matrix diagonalization method in Easyspin [6]. To simplify the analysis the nuclear Zeeman frequency was held fixed at m H = 14.7 MHz. The results in Table 3 were obtained using the Matlab functions regress and lsqcurvefit to perform fits to the data. The errors in the principal values were computed as described in [2,3]. The smallest principal value in a), obtained by linear regression and neglecting the nuclear Zeeman term, differed by 5.9 MHz from the correct value. The calculated errors in a) are therefore expected to be underestimated. The principal values and the eigenvectors obtained by a non-linear fit to the 1st order roadmaps and including the nuclear Zeeman term in Table 3b agreed exactly with those used to generate the 1st order roadmaps. The residuals were also considerably reduced, Fig. 5. The tensor data in Table 3c differed slightly from the true values, probably because the roadmaps were calculated exactly while the analysis is based on 1st order theory in this and previous works [2][3][4][5]. The slight deviations of the principal values from the true ones indicate that the high field approximation does not strictly apply in this case.
The method was originally intended for the analysis of hfc data of S = 1/2 and I=½ particles. The method is also applicable for nuclei with I = 1, provided that the nuclear quadrupole interaction is negligible and that the 1st order theory employed in [2] applies. An analysis of the X-band hfc data due to 14 N (I = 1) with principal hfc constants A N = [7, 7, 77] MHz is presented in Table 4. A fixed value was employed for the 14 N Larmor frequency, m N = 1.05 MHz.
Simulated roadmaps, accurately calculated by matrix diagonalization [6], were employed. The Matlab functions regress and lsqcurve-fit were used to obtain fits to the data. The assumed axial symmetry was approximately retained in the analysis, and the principal values were close to the correct ones, even by neglecting the nuclear Zeeman term. The small errors that do occur may in part be due to the differences in evaluating the hfc splittings, exactly in the roadmap data, to 1st order in the fits. The direction cosines of the axial component were correctly calculated. The directions of the other components were inaccurate as expected for the nearly axially symmetric tensors in Table 4.
The results support the assumption [8 pp 185-187] that the nuclear Zeeman term might be neglected in the analysis of the EPR data due to the hfc of 14 N and other nuclei with low nuclear g-factors, at least at X-band microwave frequency and lower.

The Nelson method: Analysis of EPR and ENDOR spectra of organic radicals
The method presented in [3] involved a simultaneous fit to the experimental data in three crystal planes. A FORTRAN code (Magres) was developed to obtain the g-and hfc-tensors of organic radicals by a linear fitting procedure including an analysis of the probable errors.
The tensors were assumed to be symmetric with components calculated from Eq. (11) by a standard linear regression analysis.
The predictor functions X 1-5 in Eq. (11) were expressed in terms of the polar angles specifying the direction cosines (l x l y l z ) of the magnetic field in the LAS. The parameters A 0-5 were estimated by simultaneous fitting to the experimental data in three crystal planes.
Hfc tensors were determined in a similar way in single crystal ENDOR measurements by fitting to the quantity m 2 AE À m 2 N where m N is the nuclear Zeeman frequency, m ± the ENDOR frequencies corresponding to m S = ±½. A set of angles was introduced in a final step to allow for the experimental uncertainties in the direction cosines and optimized by a non-linear procedure. The software is also applicable for obtaining the hfc-tensor from EPR data, provided that the nuclear Zeeman energy is negligible, otherwise the procedure by Fouse and Bernhard [2] might be used. The Table 3 Hfc tensor of an a-H proton [19] a) excluding, b-c) including the nuclear Zeeman energy. a) Using Matlab function regress to analyze 1st order roadmaps b) Using Matlab function lsqcurvefit to analyze 1st order roadmaps, c) Using Matlab function lsqcurvefit to analyze exact roadmaps. method was adopted in several applications involving EPR measurements of g-tensors and ENDOR studies of organic radicals with hfc due to nuclei with I = ½ [19,20,26]. Equivalent expressions were employed in the present work with n = 2 for the g-factors, hfc-splittings and ENDOR frequencies.
Based on the Hamiltonian H ZFS = S.D.S in the LSS, Eq. (13) was applied for the zero-field splitting F in an arbitrary direction (l x l y l z ), assuming that the high field approximation was valid (and the g-anisotropy negligible).
Note that n = 1 in Eq. (11) for this case and that the condition Trace(D) = 0 was not imposed.
The hfc tensor was obtained from the high frequency branch of the ENDOR roadmaps in three crystal planes. The principal hfc constants A i in Eq. (14) with errors dA i were calculated by assuming that the roadmaps corresponded to the m S = -½ state. The eigenvalues t i of the tensor T ¼   Table 5 Calculated principal values (MHz) and directions of a) an a-H and b) a weak b-H coupling using three different procedures.
The procedure was adopted from the Magres software [3] and implemented in the Nelson as well as the Schonland functions. Simulated roadmaps of an 1 H hfc tensor typical for organic radicals were analyzed in order to validate the method. The principal values calculated with (11) agreed with the assumed values A = [À5 30 40] MHz to within four decimal places; the direction cosines were also in excellent agreement with the assumed ones. The case with t i < Àm 2 N can at present not be handled by the Nelson and Schonland functions, however.
The orientation of the sample can in general not be observed during measurements and the angle h 0 in Eq. (10) can be difficult to determine precisely. We therefore assumed that the angles h 0 at B||y, z and x upon rotation about x, y and z, respectively, should be treated as parameters, employing Eq. (15) and Eq. (16) in a nonlinear fit. The corrected field direction h in Eq. (16) is expressed in terms of the uncorrected direction h 0 . Following the convention in [1], I 0 is along x, y and z, respectively, under rotations about each of the z, x and y, axes.
The equations are equivalent to Eq. (10) as verified by inserting Eq. (15) in Eq. (16) and applying trigonometric and tensor identities. The procedure was numerically validated by observing that h 0 angles deviating from the true values by a few degrees were automatically corrected by the fitting to simulated roadmaps with known parameters.
Performance: The performance was initially studied by analyzing simulated ENDOR roadmaps of the 1 H hfc due to an a-H coupling [19], and a weaker coupling [20] with the assumed tensors in Table 5. The roadmaps were calculated exactly with Easyspin [6], and by 1st order theory [22]. In the latter case the tensors computed with the Nelson function were identical to those used to generate the maps and are not presented. The tensors obtained with the Magres program [3] using the Easyspin roadmaps were included in Table 5 for comparison with those computed by the Nelson function. The tensors obtained by the Matlab functions regress, robustfit and lsqcurvefit employed in the Nelson program agreed closely; the data displayed in the table were calculated by fits with the regress function. The calculated tensors were also in good agreement with those computed by the Magres program, indicating a satisfactory performance of the Nelson function for the studied cases. The deviation of the principal values from those assumed to obtain the roadmaps is larger for the a-coupling than for the weaker coupling, indicating that the deviations are due to the inaccuracy of the 1st order theory used in the Nelson function as well as in the Magres program. One might conclude that the analysis is slightly inaccurate at X-band for the relatively large proton couplings that may occur in aliphatic radicals.
The fits to the roadmaps in Fig. 6 of the weaker coupling in Table 5 yielded a tensor in good agreement with the assumed tensor indicating that the high field approximation applied in this case.
The performance was further examined by reanalysis of the experimental g-factor data of the Stable Alanine Radical (SAR) [18]. The fitted curves (+), using the regress, robustfit and lsqcurvefit functions in Fig. 7 were slightly shifted in angle compared to experiment, indicating that the angles h 0 were not exactly positioned, as displayed more clearly on the computer screen at run time. The fit was improved by simultaneously optimizing the three h 0 angles and the six tensor elements using the lsqcurvefit function. The obtained principal values were only slightly affected by this treatment and agreed well with those obtained by the Magres program [18] and the Schonland function (Table 2), while the direction cosines were somewhat changed as seen in Table 6. Fig.6. Fits (+) to ENDOR frequency data (o) due to a weakly coupled proton by four methods. The Case number in each subplot represents data for the yz (left), zx (centre) and xy (right) planes of successively increasing crystal orientation angles (5°intervals) with respect to the external magnetic field in the three rotation planes.
An attempt was made to reanalyse the experimental ENDOR data of the a-H 8 proton of radical III in X-irradiated L-asparagine [19]. The fit was improved by simultaneously optimizing the three h 0 angles and the six tensor elements using the lsqcurvefit function.
The principal values and direction cosines in Table 7 were somewhat affected by this treatment. The roadmaps of the experimental, (o), and fitted ENDOR data, (+) shown in Fig. 9 were less suitable to elucidate the quality of the different fits because of the lower graphical resolution. It appeared, however, that the experimental and fitted data overlap most closely by employing the lsqcurvefit function with optimization of the angles h 0 in Eq. (10).
The performance was further investigated by reanalysing the experimental ENDOR data of irradiated tartrate crystals. For the decarboxylated Rochelle salt radical R2 [25] in Table 8a the residuals were significantly reduced by the simultaneous optimization Fig. 7. Fits to the g 2 -factors of the Stable Alanine Radical (SAR) [18] by four different methods. The Case number in each subplot represents data of successively increasing crystal orientation angles with respect to the external magnetic field in the yz (left), zx (centre) and xy (right) mutually orthogonal rotation planes.  Table 7 Hfc tensor of the a-H 8 proton of radical III in X-irradiated crystalline L-asparagine [19]  of the tensor elements and the angles h 0 , using the lsqcurvefit function. The tensor in Table 8b was similarly obtained by analyzing the data of radical R1(3) in irradiated ammonium tartrate [20].
The hfc-tensor in Table 8b agreed exactly with the tensor obtained by the Schonland function using the same experimental data of the R1(3) coupling. The direction cosines in the table were also in good  agreement with previous results [20]. The magnitudes of the principal values were in fair agreement with those reported previously although the signs were reversed. The absolute signs were not determined in this work, however.
Procedures to further improve the performance of the software are discussed below.

Discussion
The three Matlab functions presented in this work for the determination of the g-, hfc-and zfs-tensors were prepared mainly for the analysis of data obtained by EPR and ENDOR measurements of organic free radicals and triplet molecules trapped in single crystal matrices. The software was based on established theory [1][2][3] and validated by analysis of experimental data and roadmap data simulated from known magnetic couplings. By using the available functions for regression and error analysis the number of program lines was reduced to ca 200, i.e. by an order of magnitude in comparison with the older Fortran and Basic programs used for this purpose [2][3][4][5].
Code (I) was modelled according to Schonland's method [1], which was the first usable technique to determine the g-tensor of a paramagnetic complex in a monoclinic crystal by least squares fits to the data in three crystal planes. For crystals with an orthorhombic space group the monoclinic angle was set to b = 90°. An option to obtain the tensor by rotation about the axes of the non-cartesian monoclinic CAS was employed for the analysis of experimental data like those in Fig. 3. Slight extensions were also necessary to obtain the zfs-tensor by EPR and the hfc-tensor by ENDOR measurements. The separate analysis of each crystal plane in this method was advantageous in the special case of local axial symmetry. The distance between two SO 3 À radical anions forming a pair (S = 1), was estimated from measurements in a single plane of the line splitting attributed to the magnetic dipolar interaction between the radicals. The analysis of the EPR roadmaps due to a strongly anisotropic hfc of a nucleus with I = ½ required a non-linear treatment in the case when the nuclear Zeeman interaction is not negligible, like for the a-proton couplings in p-electron radicals, [26], see also [8 pp 175-181] and references cited therein. Code (II) was prepared following the method used by Fouse and Bernhard [2]. It might be employed in place of previously available software [4] to obtain the hfc tensor of I=½ nuclei.
Linear fitting was employed in code (III) following the method in [3] to obtain the g-, hfc -and zfs -tensors by analysing EPR and ENDOR roadmap data. A simplified procedure to correct crystal misalignment was applied in certain cases to improve the fitting to experimental data by a non-linear method. Data were supplied interactively in the three codes employing tools available in Matlab.
The error analysis presented by Fouse and Bernhard [2] was further developed in [3,5] to account for the uncertainties in the orientation of the rotation axes and of the angle of rotation about each axis. The former corrections were minimized by crystal alignment employing X-ray methods [3b,25] and were at present not accounted for. The uncertainty of the angle of rotation was, however, retained in the Fouse/Bernhard and Nelson functions by including an additional parameter h 0 in Eq. (10). The treatment resulted in better fits to the experimental data in several cases, indicating that the h 0 value was not always accurately located experimentally. The treatment is not applicable in the Schonland function, where each roadmap is analysed separately. This limitation can be partly overcome by the analysis of the maxima and minima of the roadmaps in three orthogonal crystal planes [1].
Accurate modelling is particularly important in ENDOR measurements, often with measurement errors of 0.02 MHz or lower, while the 2nd order corrections might be an order of magnitude larger under unfavourable conditions [27,Toriyama]. The computations in this and earlier works [2][3][4][5] were based on 1st order perturbation theory, rather than by the matrix diagonalization technique employed in modern software [6,7]. Detailed procedures to obtain these tensors by the latter technique seem to be unavailable in the documentation [6], while the method used here was close to that originally applied for the g-tensor [1]. On one hand only slight modifications were thus required to obtain the hfcand zfs-tensors [2][3][4][5]8], on the other the development of software based on matrix diagonalization would have been an elaborate task even by employing the tools of modern software [6]. Indeed, in principle, writing the programs using exact diagonalization is nowadays quite straightforward. However, the fact that these options do not yet exist, indicate that, in practice, it is not so simple to realize such general (albeit still with limitations, e.g. high spin systems) programs with a broad applicability as has the software discussed in the present paper. It was nevertheless necessary to estimate the loss of accuracy for some typical examples. Errors of the order 0.01 MHz were previously estimated in the analysis of ENDOR data for principal values up to ca 20 MHz [8 pp 193-194] in agreement with the data of Table 3. These errors may be neglected in comparison with other errors considered below. Larger errors of the order 0.2 MHz were obtained for the hfcsplitting of methyl radicals [27]. The shift might not be neglected in single crystal ENDOR but be insignificant in EPR due to the larger linewidth. The zfs of organic triplets are usually larger [8]. According to the formulae in [23], see also [9 pp 197-199] the observed splitting for S = 1 is correct to 2nd order, however. The items considered below might also help to improve the accuracy in the determination of the tensor values.
Modeling of ENDOR resonance frequencies by 2nd order perturbation theory Adjustment of ENDOR frequencies obtained at different magnetic field settings Table 8 Hfc tensors (MHz) in irradiated tartrate crystals. a) a-H hfc of decarboxylated Rochelle salt radical R2, a*=a [25]. b) b-H hfc of ammonium tartrate radical R1(3) [20]. Influence of g-anisotropy on the values of the hfc -and zfstensors Influence of crystal misorientation on the values of the g-, hfcand zfs -tensors Equations to treat the first case have been obtained for the ENDOR transition energies correct to 2nd order of a species with S = 1/2, I=½ taking the g-anisotropy in account [23]. The formulae in [8,Eqs. 5.3.11 and 5.4.23], adapted to yield the transition frequencies to 1st and 2nd order may serve to improve the fit.
The second issue was taken in account in the Magres program by the late William Nelson [3]. The variation of the nuclear Zeeman frequencies, depending for instance on changes of the B obs value during the recording of an ENDOR roadmap, was considered. The issue was treated in a different manner by Magoub et al in an ENDOR study of the SO 3 À anion radical interacting with a 133 Cs nucleus. Since the ENDOR frequencies were determined at different B obs values, the nuclear Zeeman frequency differed slightly. The ENDOR frequencies were therefore adjusted to a common nuclear frequency m 0 N employing Eq. (17).
A simple expression for the derivative @mAE @m N was obtained when the m + and the mfrequency branches were experimentally determined [28]. The third issue was considered in an early paper by taking in account the effect of g-anisotropy on the EPR splitting due to zfs and hfc interactions [29]. The nuclear Zeeman splitting was ignored, however. The option was therefore not implemented. The effect of g-anisotropy in ENDOR was further discussed by Atherton [8, pp 187-189] for the case when the high and low frequency branches were both determined.
Several treatments have been suggested to allow for the effects of misorientation of the crystal [1, 3, 12 pp 206-210]. This issue is of interest in order to accurately establish the electronic and geometric structure of the studied species.
Another type of error may occur in practical applications of the Schonland, Fouse/Bernhard and Nelson methods [1][2][3]. The crystal might for instance unintentionally be rotated in the wrong sense (or equivalently, so-called site splitting is occurring [31]). One can show that this ambiguity gives rise to only two distinct sets of principal values [1]. A procedure to resolve this so-called Schonland ambiguity by additional EPR measurements was proposed. A simpler method is applicable in ENDOR when the high and low frequency branches of an S = ½ species are both observed. Only one of the two tensors was in general compatible with the observed angular dependence of the two frequencies, which resolves the Schonland ambiguity [30]. The issue may be further analysed when the g-anisotropy is negligible in which case Eqs. (18) and (19) The quantities (A 2 ) xx , (A 2 ) yy , A xx , and A yy are not affected by a wrong sense of rotation (or wrong choice of site) while A xy and (A 2 ) xy would both be obtained with opposite signs. The value of (A 2 ) xy in Eq. (20), obtained by matrix multiplication differs from -(A 2 ) xy in the sign of the 3rd term, however.
The obtained tensors are therefore not identical in the case the rotation was in the wrong sense in a single plane. We refer to the literature [30] for a discussion of the applicability and limitations of the method from a numerical point of view.
The programs presented here might be of interest in applications concerned with radiation research, biophysics, biochemistry and other fields where free radicals and triplet state molecules in solid materials are mainly observed. However, a further development of the software to include an analysis based on a more general Hamiltonian like in [6,7] seems desirable, in order to establish a common standard for the determination of tensor quantities.

Summary
Three Matlab functions based on established theory [1][2][3], (I) Schonland, (II) Fouse/Bernhard, (III) Nelson have been validated by analysis of experimental data and simulated roadmap data from known magnetic couplings. By using the available functions for regression and error analysis the number of program lines was reduced to ca 200, i.e. by an order of magnitude in comparison with the older Fortran and Basic programs in Table 9. The source and executable codes of the MAGRES and ENDPAQ programs are still available, the other are either unavailable or in obsolete code, indicating that software in traditional languages may no longer be maintained, while codes developed in specialized laboratories are not easily available. The software presented in this work might therefore be of interest for the analysis of EPR and ENDOR single crystal measurements, particularly for organic paramagnetic species with S=½ and S = 1 states. Code (III) is generally applicable, codes (I) and (ii) for special cases described in the text.

Source codes
The source codes in Ref. [33] can be downloaded from Mendeley Data. To run the codes the MatLab and EasySpin [6] software must also be available. Additional information is available at https://old. liu.se/simarc/downloads?l=en.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.