3.1.1. Savitzky–Golay filter

Savitzky–Golay filtering (Savitzky & Golay, 1964) involves a convolution operation that fits a subset of points (w) from the equidistant raw data with a polynomial of degree p. The central point of the subset, which must be odd, is used to calculate the new value of the ordinate and the window moves to the next point for the next iteration. The polynomial fitting has well defined coefficients based on the size of the window (w) and the degree of the polynomial (p) (Guest, 2012). This convolution procedure is known as a moving average because the central point is estimated by weighting the surrounding points for the convolution coefficients. By selecting an appropriate window of points w, the smoothed central point of the window can be treated with convolution coefficients C_i according to

The result of this operation can be observed in Fig. 2, where the red curve is obtained using the Savitzky–Golay smoothing on the points of the black curve with the coefficients described in equation (1). An additional parameter m is included in the algorithm to establish the order of derivation of the polynomial fitting function, enabling the user to calculate the numerical derivative of order m of the smoothed spectrum while smoothing the data profile. The Savitzky–Golay filtering can be included in the software in two ways: (i) by creating tables with the convolution parameters for each combination of m, p and w; or (ii) by calculating the fitting polynomial on the points at each step of the iteration starting from zero and using ordinary least squares in the matrix form:

where a is the vector containing the coefficients C_i of the interpolation as they would be obtained from equation (1), M is a Vandermonde matrix suited for building the model of the desired p degree and y is the vector of length w containing the original intensities to be interpolated. The latter method involves a slow step such as the matrix inversion, but the ROOT framework fitting procedure is very fast and efficient, and enables polynomial coefficient refinement and the fitting of the points of the window w without using the matrix notation of equation (2) relying on the Cholesky decomposition method (Benoit, 1924) included in the MINUIT minimizer. The polynomial is stored by means of the ROOT TF1 class (https://root.cern.ch/doc/master/classTF1.html) and its coefficients are available for further analysis, such as derivation.

Figure 2 Example application of Savitzky–Golay filtering of the PXRD profile of calcite. The profile before (black line) and after (red line) application of the filtering with the following parameters: order of derivation = 0, polynomial degree = 2, window = 9.

3.1.2. Peak alignment

During in situ PXRD experiments, external stimuli such as temperature or pressure variations could result in unit-cell expansion or shrinkage, thus altering the cell parameters. However, these changes may not be the focus of investigation, which may instead be directed towards structural rearrangements due to ion adsorption or the onset of new interactions. In this context, peak shifts resulting from crystal cell variations could hide interesting data trends in the diffraction signal related to the structural changes. Transforming diffraction patterns so that these peak shifts are reduced is a challenging task that involves stretching or deforming the independent variable (2θ) axis in each profile. Two different approaches were developed by Guccione et al. (2018b) as MATLAB scripts to solve the problem.

The first approach, called peak correspondence, entails detecting the peaks in various profiles and tracing of the corresponding peaks through the profiles. The shifts among corresponding peaks are then used to estimate the optimal 2θ stretch model by means of a spline interpolation.

The second approach, `harmonic time stretch', is inspired by dynamic time warping (Rabiner et al., 1978). This approach assumes that any change in the 2θ axis can be represented by a generic function. To extend the 2θ stretch model, the function is allowed to be composed of a superposition of harmonic functions, whose number, amplitude and phase are unknown. To determine the optimal unknown values, a cost function is minimized, which is the correlation coefficient of all the profiles with a reference profile that coincides with the first uploaded file.

The two approaches have been implemented in RootProf and can be used separately or in sequence. The efficiency of the peak-alignment procedure is evaluated with a figure of merit (FOM) that measures the degree of peak misalignment. This is defined as the average of the absolute difference between the peak position and the reference position, weighted by the peak height. This evaluation is carried out for all selected peaks and profiles.

Fig. 3 depicts the effect of the peak-alignment procedure on PXRD data for the first steps of a solid-state non-isothermal reaction between tetra­cyano­quinodi­methane and fluorine which leads to co-crystal formation in the final state (Guccione et al., 2018b).

Figure 3 PXRD patterns collected during a temperature-dependent in situ experiment (a) before and (b) after application of the peak-alignment procedure consisting of peak correspondence followed by harmonic time stretch, zoomed-in on a restricted 2θ interval.

3.2.1. Semi-quantitative analysis by PCA

A fast semi-quantitative analysis that can be applied to PCA results has been implemented to directly quantify components in mixtures. It is based on the barycentric coordinate system, which is a coordinate system defined within a simplex, typically useful when describing mixtures (Mobius, 1827). Each point inside the simplex can be described by a set of q coordinates that is equal to the order of the simplex minus one. The sum of all the coordinates of a point is always equal to 1, i.e. . The barycentric system is closely linked to the Euclidean coordinate system, and points can be transformed from one system to the other relatively easily. The book Experiments with Mixtures (Cornell, 2011) describes a convenient method for projecting the points belonging to a q-dimensional simplex-centroid design in a q − 1 independent coordinate system. The inverse problem, which involves passing from a set of points represented by Cartesian coordinates to the space of a simplex, can be a relatively simple operation if the set is arranged optimally in space. However, if the group of points does not correspond to a simplex projected in Euclidean space, with the right orientation along the bi­sector of the axes, the problem can be more complex and can be addressed using PCA. The orthonormal basis of the PC space guarantees that the Euclidean metric can be used to transform simplexes that may appear among the scores or between the loadings to a barycentric system for the estimation of the relative distances between one point and the others.

Code has been developed for the case q = 3 and has been successfully applied to the results reported by Lopresti et al. (2022). The algorithm uses the q − 1 generated scores as coordinates for the points, as shown in Fig. 4. The scores are all range-scaled in [0, 1], and one vertex, the nearest one, is centred at (0; 0). The q vertices of the simplex are identified, and the q − 1 PC scores are range-scaled for the proper interval indicated by Cornell (2011) as a + b. An iterative loop then searches for the correct orientation of the simplex in space by multiplying the q − 1 scores matrix by a rotation matrix. The loop refines the rotation angle by a dichotomic search until the difference between the coordinates of the vertices reaches the stop condition, which is met when the difference between the coordinates of the vertices is lower than 10⁻¹⁰. Finally, the estimated weight fractions are calculated by multiplying the PC scores, in the form of a matrix, by the inverse of the A matrix, defined as

where I is the identity matrix and J is a (q − 1)-dimensional square matrix with all elements equal to 1. The results of this operation are reported in Fig. 4.

Figure 4 Score plot of the first two PCs for ternary mixtures composed as described by Lopresti et al. (2022). Weight fraction estimations for the three crystal phases obtained by semi-quantitative analysis are reported close to each representative point. The percentages of the total data variance explained by the first two PCs (PC1 and PC2) are reported on the axes.

3.2.2. Optimal constrained component rotation

PCA is a data-reduction method used to extract general trends in large datasets. When applied to PXRD profiles collected from in situ experiments, it allows the user to describe the whole data matrix in terms of a few PCs. The score values of these PCs capture the trends over time, while the loading values indicate the contribution of each individual peak in the profile to the particular component (Guccione et al., 2021). An example is given in Fig. 5, where a data matrix composed of 25 PXRD profiles [Fig. 5(a)] is described using three PCs, accounting for 99.9% of the total data variance. Their scores [Fig. 5(b)] define the main trends in the data, which are distinct among the different components (since they are orthogonal). The scores of the first PC (PC1, which accounts for 97.3% of the total data variance) exhibit a symmetric behaviour that mainly depends on the intensity variations of the basal peak at 2θ = 3.1°, and in fact the same peak is prominent in the PC1 loadings [Fig. 5(c)]. On the other hand, the scores of the second PC (PC2), which describes 2.3% of the total data variance, show a monotonic trend. By analysing the corresponding loadings, we can deduce that this is due to a systematic shift of the peak position during the in situ experiment. In fact, the PC2 loadings have non-zero values corresponding to the 2θ values of the main PXRD peaks, and for each peak, they exhibit a negative peak slightly shifted at lower 2θ values followed by a positive peak slightly shifted at higher 2θ values. This typical shape arises from the convolution between two peaks translating with respect to each other. The third PC (PC3), which only accounts for 0.3% of the total data variance, could be considered noise given the loading values, which have similar fluctuating contributions from all the peaks of the PXRD profiles, regardless of their height. However, the PC3 scores exhibit a well defined trend, which is approximately the square of the PC1 scores. Correlating PC1 and PC3 scores produce a nice parabola [Fig. 5(d)], although its axis is not perfectly parallel to the y axis. A theory developed by Chernyshov et al. (2011) can interpret these findings in terms of kinematic diffraction arising from the active atoms, i.e. atoms varying in phase with the stimulus received in the course of the in situ experiment, while the bulk atoms in the crystal sample remain silent. The signature for this phenomenon is a relationship between two trends in the data, where one is squared with respect to the other. According to this theory, called modulation enhanced diffraction (MED), the diffraction signal related to the squared component should only have contributions from active atoms (Caliandro et al., 2012). Following these arguments, a procedure to constrain PCA on the basis of the MED theory was developed, where two orthogonal PCs previously identified are individually rotated to maximize the FOM (FOM_scores) defined by the Pearson correlation coefficient between the scores of the second component and the square of the scores of the first (Caliandro et al., 2015). This algorithm was named optimal constrained component rotation (OCCR) and its application is illustrated in Fig. 5(e). Once the standard PCA protocol is completed, by activating the MED theory option, RootProf searches for the pair of components related by the quadratic relationship in the scores. In the present case study, PC1 and PC3 are rotated by the angles θ and ψ identified by OCCR. The new scores and loading values are shown in Figs. 5(f) and 5(g), respectively. The correlation between the rotated values of the PC1 and PC3 scores produces a parabola with the axis perfectly aligned with the y axis [Fig. 5(h)]. Notably, a dominant peak at 2θ = 3.1° appears in the rotated PC3 loadings, which is the supposed diffraction pattern related to active atoms (note that a sign ambiguity is present in the PCA loadings). At the end of the MED procedure implemented in RootProf, two figures of merit are calculated: FOM_scores and the normalized difference between the positive and the negative parts of the area underlying the second loading [denoted by FOM_loadings by Caliandro et al. (2015)]. Their values go from 0.89 and 0.08 before OCCR to 0.96 and 0.25 after OCCR, respectively (both figures of merit have a maximum value of 1, which indicates the best performance). The loadings of the rotated second component can be written in an output file and used as input for the crystal structure solution of the substructure, which is composed of the active atoms, as was demonstrated by Palin et al. (2015).

Figure 5 Steps involved in the OCCR applied to in situ PXRD measurements. The data matrix formed by the collected PXRD profiles (a) is processed by PCA, obtaining (b) the scores and (c) the loadings of the first three PCs. (d) Scatter plot of the scores of the first and third PCs, with the percentages of the total data variance explained by the two PCs reported on the axes. (e) FOM explored by the OCCR procedure, as a function of the θ and ψ rotation angles, where the initial and final rotations are indicated by red dots. (f) Scores and (g) loadings plots of the first and third PCs, recalculated by applying the best rotation found by OCCR. (h) Scatter plot of the recalculated values of the scores of the first and third PCs.

3.2.3. Kinetic analysis

A common issue encountered when processing data collected from in situ/operando experiments is identifying the reaction coordinate, i.e. a one-dimensional coordinate that signifies the advancement of the sample on a reaction pathway. One approach to resolving this issue is through PCA, where the scores of the first PC can closely approximate the reaction coordinate as it represents the main trend in the data. RootProf offers a solution for visualizing and writing scores of PCs in output files using the writescore command. The scores are sorted on the basis of the total data variance they explain. Once the user has assessed that the scores of a specific PC (usually the one with the highest eigenvalue, referred to as PC1) can be used to represent the reaction coordinate of the studied sample, such scores are stored in the output file (named ScoresPC1) and RootProf can be re-run to perform the kinetic analysis. To better visualize the steps involved in this process, we show in Fig. 6 the most relevant RootProf plots resulting from a kinetic analysis applied to PXRD profiles collected by a laboratory diffractometer equipped with a system for measurements at a controlled temperature on a Pd-based vapochromic compound [see Belviso et al. (2021) for further details, referring to the compound Pd(ppy)]. The set of PXRD patterns collected at different temperatures [Fig. 6(a)] is subjected to PCA, which determines the PC1 scores. In this case, the first PC explains 96% of the total data variance [Fig. 6(b)]. These are rescaled in the range [0, 1] and plotted against sample temperature [Fig. 6(c)]. The masterplot technique (Goiss, 1963) is used to choose the best kinetic model that most accurately describes the experimental data [Fig. 6(d)]. A least-squares procedure [Fig. 6(e)] is applied to fine-tune the best kinetic model, which in this case is A3 (i.e. the Avrami–Eroféev model), and determine the order of reaction. Finally, the kinetic parameters are estimated from the plot reported by Coats & Redfern (1964) [Fig. 6(f)]. For more information on the kinetic analysis, refer to Guccione et al. (2018a). The graphics in Figs. 6(c)–6(f) are generated by the second run of RootProf using the command file reported in Section S5 of the supporting information.

Figure 6 Steps involved in the kinetic analysis of in situ/operando measurements. (a) The data matrix formed by the collected PXRD profiles is processed by PCA, (b) obtaining the PC1 scores. (c) Scores are rescaled to range between 0 and 1 and are plotted as a function of the temperature on the sample. (d) A masterplot is generated to compare experimental data with a set of kinetic models and (e) the reaction order parameter for the chosen model is determined by fitting with experimental data. (f) The Coats & Redfern (1964) plot allows us to determine the activation energy and frequency factor parameters by fitting with experimental data.

3.3.1. Multithreading

Parallelization refers to the process of breaking down software operations into smaller tasks and distributing them across multiple processors or cores. This can be achieved either implicitly through multithreading or explicitly by dividing the tasks between processors. A recent study by Piparo et al. (2017) examined the differences between these two approaches to parallelization within the context of the ROOT framework and showed how the two methods can be complementary. The main goal of parallelization is to improve the performance of computer systems and reduce the execution time of an algorithm. The RootProf parallelization implementation currently utilizes implicit parallelism, and it is not yet optimized for computer clusters. However, the computational speed has improved significantly on the host computer architecture. This is evident in the quantitative analysis module which achieved a performance boost of almost a tenth of the previously required time on a machine powered by an AMD Ryzen 5 3600X processor with 6 cores, where each core can handle 12 threads.

3.3.2. Standard addition method

Quantitative analysis of the weight amount of an individual crystal phase in a complex mixture of unknown composition can be challenging. In the work of Zappi et al. (2019), the quantification of the active pharmaceutical compound (API) paracetamol (polymorphic form I) in a solid formulation of the commercial drug Tachifludec was considered as a case study. In the first instance, the commercial sample was measured by PXRD in three replicates and the known crystal phase was analysed using the MultiFit procedure of RootProf, as shown in Fig. 7(a) (samples No. 0, 1, 2 and inset on the left). The weight fraction averaged on the three replicates was 0.137 ± 0.012, which provided a preliminary estimate of the amount of paracetamol present in the commercial drug. To obtain a more accurate measurement, the standard addition method (SAM) was employed, whereby Tachifludec samples with known amounts of the API were measured. MultiFit was applied to these additional profiles, and the weight fraction estimates are shown in Fig. 7(a) (samples 3–11, inset on the right). This information was then used in combination with the known additions to produce the calibration plot shown in Fig. 7(b). The SAM determination of the API concentration was found to be 0.141 ± 0.010 by calculating the intercept of the regression line on the y axis, using samples 3–11. Alternatively, the weighted mean between the intercept of the regression line on the y axis and the absolute value of its intercept on the x axis was considered, which was found to be 0.152 ± 0.009. Notably, the nominal concentration of the API in Tachifludec is 0.145. RootProf performs these calculations by providing the known added concentration of the standard as input, as reported in Section S7 of the supporting information.

Figure 7 Quantitative analysis to assess the amount of paracetamol present in Tachifludec. (a) Weight fractions estimated from powder samples of Tachifludec, where paracetamol has been added at concentrations of 0, 15, 25, 35%(w/w) with three replicates per sample. A linear fit has been carried out on the samples without standard added (samples 0, 1, 2). Insets show the fit of the pure paracetamol profile (red line) on the first (left) and last (right) samples (black line). (b) Calibration plot obtained by the SAM procedure, where a linear fit has been carried out on the samples with the standard added (samples 3–11). Errors are calculated from the whole-profile fitting procedure on single mixtures, and regression lines are superimposed in red.

3.3.3. Combined fitting of different datasets

Quantitative estimates by PXRD may be difficult to obtain when reference profiles have limited defining features. In fact, achieving a reliable quantitative phase analysis is quite challenging for poorly crystalline, nanosized or quasi-amorphous mixtures. This limitation has a notable impact on the pharmaceutical field, especially given that amorphous compounds can potentially offer higher bioavailability or may be the only option when crystal products cannot be obtained (Kavanagh et al., 2012). However, the PDF approach provides a convenient means of investigating these types of compounds. By interpreting total scattering in direct space and defining the probability G(r) of finding any two atoms at a given interatomic distance r, the PDF allows for more accurate qualitative analysis (Egami & Billinge, 2012). The MultiFit approach of RootProf has been adapted to improve the quantitative analysis of low-crystallinity compounds, by combining the results obtained in direct and reciprocal space. The approach, as shown in Fig. 8, involves collecting PXRD patterns on nine binary mixtures containing microcrystalline cellulose and the Eudragit L100 polymer (Mahdi, 2015) using a PDF setup that required the sample to be placed very close to the detector to acquire data at higher momentum transfer. A linear combination of the two pure-phase profiles was used to fit the data separately in reciprocal space [Figs. 8(a) and 8(b)] and in direct space [Figs. 8(c) and 8(d)]. By comparing Figs. 8(b) and 8(d) we can observe that Eudragit L100, which represents the quasi-amorphous phase (phase 1), was better quantified in reciprocal space, whereas MCC, which has higher crystallinity, was better quantified in direct space. Quantitative estimations in the direct and reciprocal spaces were then averaged, resulting in the improved estimates shown in Fig. 8(e), where the combined contribution from the dual spaces provided equally precise estimations for both phases. The Kullback–Leibner distance (Kullback & Leibler, 1951) was used to quantify the improvement obtained, as reported in Figs. 8(b), 8(d) and 8(e).

Figure 8 Quantitative analysis of a poorly crystalline pharmaceutical mixture, processing PXRD data in reciprocal and direct space. Least-square fits in (a) reciprocal and (c) direct space for sample number 8, are shown, with observed and calculated profiles shown in black and red, respectively. Estimated weight fractions of Eudragit L100 (black, Ph 1) and MCC (red, Ph 2) determined in (b) reciprocal and (d) direct space, and (e) obtained as the mean value between the two spaces. Error bars indicate fitting errors in (b) and (d), and propagated errors in (e). The sum of the weight fractions for each sample (in green) and the values of the Kullback–Leibner distance (AKLD) are reported.

A step forward in combining direct and reciprocal space datasets was attempted by implementing a procedure for combined fitting, but worse results were obtained, possibly due to the difficulty in properly weighting the profiles across the two spaces.

Notably, performing pre-processing of data could lead to fundamental improvement in quantitative analysis, and this could differ for the two datasets. With this aim, the RootProf procedures for supervised quantitative analysis, employing all possible pre-processing options sorted by a proper FOM, can be used to select the best pre-processing option.