1 Introduction

Infrared and terahertz (THz) spectroscopies are incredibly powerful analytical techniques with many applications across the physical and life sciences. Whilst the origin of many of the spectral features in an infrared spectrum can be easily interpreted with a knowledge about characteristic vibrational frequencies of functional groups, to identify the origin of all peaks (particularly below 1000 cm− 1) and understand the subtleties in peak shape and position theoretical support is essential. There are now a number of density functional theory (DFT) based packages designed for both small molecules [1, 2] and periodic solids [3,4,5,6,7] capable of calculating vibrational frequencies and infrared intensities that can be used to interpret complex experimental spectra. There are also post-processing tools such as PDielec [8] which take into account effective medium approximations [9], the attenuated total reflection (ATR) effect [10] and Mie scattering [11] to aid in the interpretation of complex experimental spectra.

Whilst these packages are now readily accessible, the calculation of spectra of complex systems, that correlate well with experiment can still be tricky, particularly at frequencies below 200 cm− 1 [12,13,14]. Choice of basis set or pseudo-potential [15, 16], density functional [17], DFT package [16] and convergence criteria can all have a dramatic affect on any calculated spectral parameter. This accessibility has a downside and can lead to unsuitable calculations being compared with experiment, for instance it is still commonplace to compare single molecule DFT calculations to THz measurements of crystalline material [18,19,20,21,22,23] which has often been shown to be unsuitable [24].

In this paper we will discuss the potential pitfalls in the calculation of the infrared spectra of crystalline materials using a range of solid state DFT packages. We will look at the reliability of such calculations with respect to method, basis set and pseudo-potential. We will also compare a number of van der Waals’ dispersive corrections which have been shown to be particularly important in the calculation of THz spectra [25, 26].

The crystal structure of the material, sodium peroxodisulfate, is available from the Crystallographic Open Database [27] with code number 2208366 and was determined from single crystal X-ray measurements at 150 K [28]. Sodium peroxodisulfate forms an interesting monoclinic, P1¯, crystal structure with half a formula unit in the asymmetric cell and one formula unit in the full cell. Each sodium ion is in a distorted octahedral environment of oxygen atoms. The two SO4 groups are linked together by a bridging O–O bond with an experimental bond length of 1.479 Å, as can be seen in Fig. 1. This is an unusual functional group which will test the suitability of any oxygen basis set or pseudo-potential. Sulphur can also be a problematic element in THz spectroscopic calculations [29] again providing a challenge for the choice of basis set or pseudo-potential. Finally, the sodium ions are bound in the lattice by ionic bonding, causing considerable charge separation which means it is likely that the infrared spectrum of the material is influenced by the crystal shape. Previous work [8, 10, 30] has shown that for ionic materials it is necessary to take into account the interaction between the crystal phonon modes and the electric field of the infrared radiation as large shifts in the absorption peaks can occur.

Fig. 1
figure 1

The unit cell of sodium peroxodisulfate—Na (purple), S (yellow), O—axis labelled a, b and c with the origin labelled o

Full size image

2 Methods

The following density functional theory (DFT) packages were used to calculate the optimized geometry and unit cell of Na2(SO4)2; Abinit [5], CASTEP [3], Crystal17 [31], VASP [4] and Quantum Espresso [7] (referred subsequently to as QE). For all, except Crystal, pseudo-potentials were used to represent the core electrons. Crystal14 [6] was used for some of the preliminary calculations reported in the Supplementary Information (SI). The Perdew-Burke-Ernzerhof (PBE) functional [32] was used for all calculations. The pseudo-potentials, basis sets and the associated packages are summarized in Table 1 along with the atomic configuration of the active electrons in the calculation.

Table 1 DFT packages, pseudo-potentials, basis sets and atomic electronic configurations
Full size table

The ONCVPSP pseudo-potentials [33] and FHI pseudo-potentials [34] were obtained from the Abinit website [35]. The CASTEP 19.1 norm-conserving pseudo-potentials were taken from the ‘on-the-fly’ pseudo-potentials built using the NCP19 keyword as input to the ‘SPECIES_POT’ directive. In some of the early calculations that have been included for completeness, CASTEP 17.2 with NCP17 pseudo-potentials had been used and this is indicated where necessary in the SI. The QE pseudo-potentials are Ultra Soft Pseudo-potentials (USPs), taken from the accurate set of Standard Solid State Pseudo-potentials (SSSP) [16, 36] and the VASP pseudo-potentials are Projector Augmented Wave (PAW) pseudo-potentials [37] distributed with VASP 5.4.4.

Crystal is an all electron code using atom centred contractions of Gaussian functions to represent the electronic wavefunction. Two basis sets were used; a triple zeta valence plus polarization (TZVP) basis [38] and a larger DEF2 basis set, based on the def2-TZVP molecular basis of Weigend and Ahlrichs [39]. The TZVP basis was taken from the library of basis sets on the Crystal website [40]. The DEF2 basis set was obtained from the Basis Set Exchange web site [41]. All f-functions were removed from the basis set and any Gaussian functions with an exponent less than 0.1 were removed from the sodium basis set in order to prevent linear dependencies in the calculation. Details of the number of uncontracted Gaussian functions and their contractions are given in the SI. During the course of the work reported here Crystal17 became available along with additional options for the dispersion correction. Crystal14 and Crystal17 gave almost identical results for calculations using the TZVP basis. It was decided therefore to report results using Crystal17 for all calculations reported in the main text. Calculations regarding the optimization of the dispersion correction parameters were performed with Crystal14 and they are reported in the SI.

Alongside the choice of pseudo-potential there are also a number of possible choices for a suitable description of the dispersive interaction in solid state DFT codes. Table 2 shows the dispersive correction used with each package;

Table 2 DFT packages and the dispersion corrections used
Full size table

The Grimme D2 method [42] adds a semiempirical correction to the energy of the form;

$$ E_{\text{disp}} = -\frac{1} {2}\sum\limits_{i\neq j}^{} \frac{C_{6ij}} {r^{6}_{ij}}f(r_{ij}) $$
(1)

where the damping factor is;

$$ f(r_{ij}) = \frac{S_{6}}{1+e^{-d({r_{ij}}/{S_{R} R_{0ij}} - 1)}} $$
(2)

The S6 and SR parameters are usually taken from the literature after fitting to reference calculations and are specific to the functional being used; for the PBE functional and the GD2 method values of 0.75 and 1.0 respectively are used. It has been common to vary their values so as to minimize the deviation between the calculated and experimental crystal structures. In some cases this is achieved by minimizing the error in the calculated volume [43, 44], in others the distortion of the cell is minimized [45].

Becke and Johnson [46] introduced an alternative damping function which leads to a constant contribution to the dispersion energy for spatially close pairs of atoms. The Grimme D3 method (GD3) [47] is a revision of the GD2 method which includes 3-body terms, allows for some geometry dependent information to be taken into account and includes higher multipole contributions. This method can also use the Becke and Johnson damping scheme (GD3-BJ).

The method of Tkatchenko-Scheffler (TS) [48] is formally equivalent to the GD2 method but the parameters are functions of the charge density. For this dispersion correction the S6 parameter is set to 1.0 and the parameter SR is used for fitting (0.94 is the default value for the PBE functional).

2.1 k-point Integration

All calculations used a k-point integration scheme based on the Monkhort-Pack method [49] with 7, 6 and 5 points in each reciprocal lattice direction respectively. Based on optimizations of the unit cell and the molecular geometry using VASP with an energy cutoff of 560 eV, increasing the k-point grid density by a factor of 2 showed that the energy of the optimized unit cell changed by less than 0.001 eV and the change in the calculated volumes of the optimized unit cells was less than 0.01 Å3. Further details are given in the SI. Calculations of the phonon spectrum at the gamma point showed a difference of less than 0.05 cm− 1 in the frequencies of the lowest nine non-zero modes. Although no explicit k-point convergence testing was performed with the other methods, it is expected that for this insulator, the choice of k-point sampling density suggested by VASP will be equally accurate for all methods.

2.2 Plane-Wave Energy Cutoffs

Using the above k-point integration grid, optimizations of the atomic positions were performed at the, fixed, experimental unit cell dimensions. From the atom-only optimized structures the internal pressure of the cells was calculated and the phonon modes calculated at the gamma point. As reported in the SI, these calculations were carried out using a number of plane-wave energy cutoffs. Table 3 shows the chosen cutoff energy for each program and pseudo-potential along with the absolute difference between the calculated pressure, unit cell energy and frequencies for the chosen cutoff and the largest cutoff used. The root mean squared shifts in frequency were calculated by taking the frequencies of all the optical phonon modes of the calculation with the largest cutoff as a reference.

Table 3 Chosen energy cutoffs and associated absolute errors in energy, pressure and root mean squared shift (RMSS) in frequency
Full size table

The results for Abinit/FHI have the largest change in the mean squared frequencies. However, this was owing to the relatively poor translational invariance of the calculation; before the translational modes were projected out, the 3 lowest frequencies were far from zero and this caused some contamination of the two lowest optical modes. The wide range of cutoff energies reflects the type of pseudo-potentials being used by the packages. Those using USP or PAW pseudo-potentials are much softer than the norm-conserving pseudo-potentials and therefore need a smaller number of plane-waves in the basis set for a similar accuracy.

Using the plane-wave cutoffs shown in Table 3 and starting with the experimental crystal structure [27] the unit cell dimensions and the atomic positions were optimized, maintaining the space group symmetry.

2.3 Calculation of Infrared Absorption Spectra

The PDielec package [8] was used to calculate the absorption spectrum from the Born charges and dynamical matrix calculated by each of the DFT packages. The dynamical matrix was calculated by first optimizing, within the constraints of the space group symmetry, the unit cell and atom positions. This was followed by the calculation the dynamical matrix without any requirement for translational invariance. The convergence settings used by each program are reported in the SI. Translational invariance was imposed using the PDielec package which applies projection operators in real-space to project out the three translationally invariant modes of motion from the dynamical matrix. By using PDielec to perform the projection and calculation of the phonon normal modes a consistent set of atomic weights were used for all of the DFT packages.

For the calculation of the effective permittivity PDielec assumes that the sodium peroxodisulfate is a powdered crystal dispersed in a supporting matrix and that the composite material has an effective complex permittivity, calculated using an effective medium theory such as Maxwell-Garnett or Bruggeman [9]. The effective permittivity is calculated from the calculated permittivity of the sodium peroxodisulfate, the shape of the crystal and the permittivity of the supporting matrix. The effect of crystal shape on the absorption can be studied by comparing the effective medium theory spectra with that calculated at low concentrations using the Averaged Permittivity (AP) method [8], which shows absorption at the transverse optical (TO) phonon frequencies. In addition, using PDielec, it is possible to incorporate the effects of scattering during the absorption process using the Mie method for spherical particles [11]. This allows the program to describe the effect of air inclusions in the sample and to describe what effect larger particles of Na2(SO4)2 may have on the spectrum.

The attenuated total reflection (ATR) technique is commonly used to record infrared spectra. An effective medium theory calculation of the effective permittivity for a high volume fraction of sodium peroxodisulfate embedded in air enables the calculation of the reflectance spectrum in an ATR configuration. The method used in PDielec to calculate the ATR reflectance is similar to that used by others [10, 50]. PDielec solves Fresnel’s equations [51] for 45° incident radiation on a slab of non-absorbing, high permittivity material (such as diamond) supporting a layer of this effective medium.

2.4 Comparison of Calculated Infrared Absorption Spectra

In order to compare the spectra calculated by the various packages, a spectrum was calculated using PDielec assuming a Maxwell-Garnett effective medium model of 10% by volume of small particles of sodium peroxodisulfate suspended in a Polytetrafluoroethylene (PTFE) matrix. Each spectrum was calculated with a frequency resolution of 0.1 cm− 1 and the width of each absorption was taken to be 5 cm− 1. In order to calculate a normalized cross-correlation between the calculated spectra, each spectrum, A, was normalized.

$$ A_{\text{norm}}(i) = \frac{A(i) - \bar{A}}{\sqrt{n}\sigma(A)} $$
(3)

where n is the number of data points in the spectrum, \(\bar {A}\) is the mean value of the spectrum and σ(A) is its standard deviation. The normalized cross-correlation coefficient can take values from − 1 to + 1, a value of 0 indicates no correlation between the spectra.

The maximum cross-correlation between two spectra is calculated at a given frequency shift. The value of the maximum correlation coefficient and its ‘lag’ is used to calculate the similarity between the calculated spectra.

Using either the full cross-correlation matrix or the matrix of frequency shifts between all pairs of calculated spectra, a heat-map was calculated along with a clustering of the calculations according to their similarity. These calculations were performed using the gapmap package [52] in R [53].

3 Experimental

Sodium peroxodisulfate (99%) was bought from Sigma-Aldrich and ground using a Specamill stainless steel ball mill. The powder was mixed with a non-absorbing matrix material (for infrared measurements KBr was used and for THz measurements PTFE was used) and pressed using 7 t of force into pellets approximately 500 μm in thickness supported by a surrounding copper ring. A Nicolet iS5 FTIR was used for the transmission infrared measurements and 32 scans recorded for both background and sample at a 1 cm− 1 frequency resolution.

For ATR infrared measurements a Bruker alpha platinum ATR instrument was used. A sample of sodium peroxodisulfate was placed onto the diamond ATR crystal, clamped in place, and an average of 32 scans with a spectral resolution of 2 cm− 1 was recorded. THz spectra were recorded on a home-built THz time-domain spectrometer (THz-TDS) previously described elsewhere [54]. In brief, spectral measurements were performed using a dry-air purged broadband THz-TDS using a mode-locked Ti:sapphire laser (Vitara, Coherent) which was used to produce a train of near-infrared pulses, each of duration 20 fs, centred at 800 nm at a repetition rate of 80 MHz. The beam was then focused onto a low-temperature-grown gallium arsenide (LT-GaAs) on quartz [54] photoconductive switch with a large-area slot electrode design which was 200 μm wide and 4 mm long. The emitter was biased at 350 V using a 7 kHz modulation frequency with a 50% duty cycle to enable lock-in detection. The THz radiation emitted from the photoconductive switch was collected and collimated from the side of the emitter excited by the laser and focused onto the sample pellet by a set of off-axis parabolic mirrors. The THz radiation transmitted through the sample was then recollected and focused with a second pair of mirrors onto a second LT-GaAs-on-Quartz device used as a photoconductive detector. The current generated in the photoconductive switch was amplified using a transimpedance amplifier with a gain of 1 × 108 Ω with a time-delayed probe beam (100 mW) split off from the original near-infrared laser pulse train used for detection. Spectra are an average of 60 scans recorded with a frequency resolution of 0.8 cm− 1. Low temperature measurements were performed by mounting the sample pellet, within a copper ring for good thermal contact, onto a coldfinger of a continuous-flow helium cryostat (MicrostatHe, Oxford Instruments) equipped with polymethylpentene (TPX) windows.

4 Results and Discussion

For comparison purposes, the experimental unit cell dimensions [28] are reported in Table 4 along with the length of the O–O bond in the crystal.

Table 4 Experimental unit cell dimensions of Na2(SO4)2
Full size table

4.1 Geometry Optimization

Table 5 shows the percentage errors in the optimized unit cell parameters for the methods without a dispersion correction. The calculated cell dimensions are provided in the SI. The calculated volumes from these optimizations are systematically larger than the experimental cell volume by more than 6%. The calculated O–O bond is also too large by more than 1.5%. This systematic error in volume is expected because the uncorrected DFT methods do not include any electron correlation, and this tends to increase the unit cell volume. There are two issues to be considered when comparing calculated and experimental unit cell dimensions. The first issue is that the experimental unit cell dimensions were determined at 150 K. For the purpose of comparison with calculation determination at a lower temperature would be better, but based on the expansions coefficients of mirabilite (Na2SO4(D2O)10) [55] the unit cell volume of sodium peroxodisulfate would only be 0.36% smaller at 0 K and each unit cell dimension would be only 0.12% smaller. The second issue is that in the calculation no account has been taken of the zero-point motion taking place in the crystal. Zero-point motion will tend to increase the volume of the crystal, values of up to 3% have been reported [56]. Whilst such calculations are now feasible, this effect is often ignored. Indeed after thermal effects are accounted for the PBE/GD3 dispersion correction already seems to over-estimate the cell volume by about 1% [57]. This is in agreement with the dispersion corrected results shown in Table 6, where the volume is generally calculated to be too large.

Table 5 Calculated percentage errors in unit cell dimensions using no dispersion correction (negative numbers indicate that the calculated value is smaller than experiment)
Full size table
Table 6 Calculated percentage error in unit cell dimensions using dispersion corrections (negative numbers indicate that the calculated value is smaller than experiment)
Full size table

A few calculations were performed, using Crystal/DEF2/GD2, CASTEP/NCP17/GD2, and CASTEP/NCP19/TS where the S6 parameter of the Grimme dispersion correction or the SR parameter of the TS correction was optimized to improve the agreement between the experimental and calculated unit cell dimensions. Details of the calculation of the optimized parameters can be found in the SI. In these cases the fact that optimized dispersion correction parameters have been used rather than the default values is indicated by ‘-v’ or ‘-r’ in the method label. A ‘-v’ indicates the parameter was determined so as to reproduce the experimental volume, a ‘-r’ indicates that the root mean squared fractional deviations (RMSFD) of the calculated unit cell dimensions and angles from the experimental values were minimized. Parameter optimization for Crystal/DEF2/GD2 was performed using Crystal14, but all phonon calculations were performed with Crystal17.

In the case of Crystal/DEF2/GD2-v, the resulting value of S6 (0.92) was larger than the default value of 0.75 and resulted in a distortion of the unit cell, as evidenced by an RMSFD of 2.8%, which is greater than that of the non-dispersion corrected unit cell (1.8%). Optimizing the parameter to minimize RMSFD, Crystal/DEF2/GD2-r, gave an S6 value of 0.5 (RMSFD 1.0%), but the calculated cell volume is 4% larger than experiment. However, the distortion is lower than that of the non-dispersion corrected unit cell.

Attempts were made to reproduce the experimental volume by optimizing the S6 parameter for CASTEP/NCP17/GD2. The results are reported in the SI and indicate a sudden change in the packing of the cell around the value needed to give the experimental volume. The value of S6 which minimized the RMSFD was the same as the default value (0.75). Although these calculations used the NCP17 rather than the NCP19 pseudo-potentials, on the basis of these results it was decided to include only CASTEP/NCP19/GD2-r and not the GD2-v results.

For CASTEP/NCP19/TS the value of SR which gave the experimental volume was 0.925. As volume and RMSFD were behaving similarly with respect to changes in SR (see SI), only the default CASTEP/NCP19/TS (SR= 0.94) and CASTEP/NCP19/TS-v results are reported.

All methods which include a dispersion correction lead to a reduction in the calculated volume of the cell. Of those methods where no parameter optimization was performed, VASP/PAW/TS and Abinit/FHI/GD2 are closest to the experimental volume. VASP/PAW/GD3-BJ and Abinit/ONCVPSP/GD2 have the lowest RMSFD (0.5%), which is lower than that achieved by Crystal/DEF/GD2-r, even with minimization of RMSFD.

The process of optimizing the S6 parameter by minimizing the error in the volume of the cell can lead to some unintended consequences. This is most clearly shown by the Crystal/DEF/GD2-v results. Although the volume calculated by this calculation agrees with the experiment, the percentage errors in the cell dimensions are much larger than that of the non-optimized case.

The O–O bond length is calculated to be larger than the experimental value for all calculations. The addition of dispersion corrections does not significantly change the value. The Crystal/TZVP and its dispersion corrected methods calculate the longest O–O bonds, indicating that there may be a problem with the oxygen basis set used for this calculation.

4.2 Translational Invariance

Accurate calculations of the phonon modes at the gamma point of the crystal should give three acoustic modes which have zero frequency reflecting the translational invariance of the crystal. If the calculation has not converged sufficiently in sampling the Brillouin zone or in basis set, or there is some underlying grid (e.g. FFT grid) being used within the calculation that is insufficiently fine, then these modes may have non-zero values. In some cases an imaginary value indicates that the cell is unstable with respect to atomic motion. No such instabilities were encountered for the calculations reported here. Any (small) imaginary frequencies were as a result of losing translational invariance. In the case of CASTEP a FINE_GRID_SCALE parameter of 6 was chosen to minimize the deviations of the acoustic mode frequencies from zero as described in the SI. Table 7 shows the root mean squared error (RMSE) in frequencies of the acoustic modes at the gamma point, as calculated by each package, without imposing translational invariance. The squared error is averaged over all the methods used by each package. The error in the acoustic mode frequencies depended principally on the package used and was independent of the method. In addition the Table shows the root mean squared shift (RMSS) of the optical frequencies for each calculation as a result of imposing translational invariance on the dynamical matrix, by projecting out the translational degrees of freedom from the dynamical matrix using PDielec and recalculating the phonon frequencies.

Table 7 The root mean squared error (RMSE) in the acoustic mode frequencies and the root mean squared shift (RMSS) in the optical frequencies on projection of the translational degrees of freedom
Full size table

Abinit shows large deviations from zero in the unprojected acoustic mode frequencies and also shows a significant effect on the optical frequencies as a result of projecting out the translational modes from the dynamical matrix. The shift in frequency as a result of projection is largest for the lower frequency modes. CASTEP, Crystal and VASP have only small deviations from zero in the unprojected acoustic mode frequencies and there is only a small shift in the optical mode frequencies after projection. The RMSS for CASTEP is 5 times lower than any other the method, which can be attributed to optimization of the FINE_GRID_SCALE parameter as reported in the SI.

The QE results show a significant deviation from zero in the unprojected acoustic mode frequencies, but this does not lead to significant shifts in the optical mode frequencies upon projection Table 7.

4.3 Calculated Frequencies and Intensities

Figure 2 shows a comparison of the non-dispersion corrected calculations of frequencies and intensities over the full frequency range. The spectrum of Na2(SO4)2 falls into three distinct ranges; low (below 300cm− 1), intermediate (300–750cm− 1) and high (above 750 cm− 1) The low frequency region up to 300 cm− 1 is complex as can be seen in Fig. 3 but all calculations predict six absorptions of varying intensity and position in this region. The calculated frequencies and intensities for all calculations are given in the SI.

Fig. 2
figure 2

Non-dispersion corrected frequencies and intensities

Full size image
Fig. 3
figure 3

Non-dispersion corrected frequencies and intensities—low frequencies

Full size image

4.4 Phonon Mode Analysis

The make-up of each phonon mode in terms of either the internal/external contributions or in terms of the contributions from particular groups of atoms can be determined from their percentage kinetic energy contribution to that mode. The approach adopted in PDielec for this analysis is described in previous work [13].

The internal contributions can be regarded as molecular vibrations and the external contributions as whole molecule translatory or rotatory motion. There are two obvious ‘molecular’ groupings for this crystal. In one the S2O8 moiety can be treated as a single unit or it can be treated as two independent SO4 units. The sodium atoms are treated as being able to contribute to external translatory modes. Figures 4 and 5 show the results of the analysis of the Crystal/DEF2 calculation using SO4 and Na molecular groups. Analysis of the results from all the other packages, pseudo-potentials and disperion corrections gives very similar results to those presented here for Crystal/DEF2. Figure 4 shows the contribution of external (translatory and rotatory) and internal (vibrational) contributions. Figure 5 shows the break-down into ‘molecular’ components. The frequencies of the individual normal modes for this particular calculation are given in the supporting information, as are the related results where the S2O8 unit is treated as the molecular unit.

Fig. 4
figure 4

Mode analysis of Crystal/DEF2 phonon modes, internal and external contributions

Full size image
Fig. 5
figure 5

Analysis of Crystal/DEF2 phonon modes, molecular contributions

Full size image

The figures show the three acoustic modes (modes 0 to 2) to be purely translational modes, as expected. All modes up to but not including mode 17 are external modes dominated by rotatory and translatory motion of the molecular groups. Mode 17 itself is a mixture of internal and external character. The external modes include contributions from the sodium atom translatory motion, as can be seen from Fig. 5. Above 321.9 cm− 1 (the frequency of mode 17) modes are all dominated by internal vibrations of the SO4 group. Inspection of the figures in the SI calculated assuming an S2O8 group show that the lowest optical mode (mode 3) has a large contribution from the vibration of the S2O8 group originating from movement of the O–O bond. Modes 28 to 31 show significant rotatory motion of the SO4 groups is associated with stretching of the O–O bond.

4.5 Infrared Spectra Determined Using a Maxwell-Garnett Effective Medium Theory

Unless otherwise stated the infrared spectra were calculated using PDielec [8] from the normal modes and Born charges using a Maxwell-Garnett effective medium theory model for 10% by volume of small spherical crystallites embedded in a PTFE matrix support. The Lorentzian line width for each absorption peak in the calculation was taken to be 5 cm− 1.

4.5.1 No Dispersion Correction

A comparison of the full frequency range for calculations involving no dispersion correction can be seen in Fig. 6. Because the spectrum is dominated by the high frequency range, Figs. 7, 8 and 9 show the same spectrum but over the high, intermediate and low frequency ranges respectively. If the calculations were fully converged in all aspects including basis set, k-point integration and grid representation of the charge and wavefunction, then it should be expected that the spectra should agree with one another. In the high frequency region (Fig. 7) this is clearly not the case. Although there is some qualitative agreement there are clear differences between the calculations. The absorptions in the region from 900 to 1050 cm− 1 are the bending modes of the SO4 units. The Crystal/TZVP results are more than 40 cm− 1 lower in frequency than any other calculation, including Crystal/DEF2. The Crystal/DEF2, QE/SSSP, Abinit/FHI and Abinit/ONCVPSP are higher in frequency than other calculations and the VASP/PAW and CASTEP/NCP19 calculations predict absorption at around 1000 cm− 1.

Fig. 6
figure 6

Non-dispersion corrected spectra using Maxwell-Garnett effective medium theory

Full size image
Fig. 7
figure 7

Non-dispersion corrected spectra using Maxwell-Garnett effective medium theory—high frequency range

Full size image
Fig. 8
figure 8

Non-dispersion corrected spectra using Maxwell-Garnett effective medium theory—intermediate frequency range

Full size image
Fig. 9
figure 9

Non-dispersion corrected spectra using Maxwell-Garnett effective medium theory—low frequency range

Full size image

The S–O stretch region near 1200 cm− 1 is more complex but shows a similar pattern, with the Crystal/TZVP calculation showing absorption at significantly lower frequencies and the QE/SSSP, Abinit/FHI and Abinit/ONCVPSP calculations showing absorption at slightly higher frequencies. Crystal/TZVP was previously shown to predict a long O–O bond, indicating a problem with this particular basis set. The deficiency in the basis set is manifesting itself with a poor prediction of the O–S–O bending and S–O stretching frequencies.

In the intermediate frequency region (Fig. 8) the pattern is not quite the same. All calculations except CASTEP/NCP19 and Crystal/DEF2 agree that there is an absorption around 680 cm− 1. The CASTEP and Crystal/DEF2 results are shifted to slightly lower frequency. In the region from 400 to 600 cm− 1 the Crystal/TZVP calculations shows absorption peaks at 480 and 510 cm− 1 which are lower in frequency than all the other calculations by about 50 cm− 1. In this frequency region all the plane-wave calculations are in agreement with each other. The Crystal/DEF2 results are similar to the plane-wave calculations but there are two distinct peaks at 520 and 530 cm− 1, instead of a single absorption around 530 cm− 1.

The low frequency regime shown in Fig. 9 is harder to unravel. There appears to be general agreement that there is significant absorption around 200 cm− 1 which comes from two strong absorptions. There is little agreement as to where the lowest frequency absorption occurs, although all methods predict some absorption below 100 cm− 1. Abinit/FHI predicts the lowest frequency absorption just below 70 cm− 1. Between 100 and 250 cm− 1 there are 3 strong peaks, the middle peak of which is the most intense and whose frequencies can shift by up to 20 cm− 1 depending on the package being used.

4.5.2 Dispersion Corrected Spectra

The effect of including a dispersion correction in the calculation of the unit cell and the phonon modes can be best seen by comparing the results of the VASP calculations with different dispersion corrections. Figure 10 shows the predicted spectra over the full frequency range. The intermediate and high frequency show small changes on including a dispersion correction and can be found in the SI. The uncorrected spectrum shows a single, slightly more intense absorption at about 1200 cm− 1, which comes from two transitions with similar frequencies. The dispersion corrected methods show two peaks at slightly higher frequencies.

Fig. 10
figure 10

IR spectra from VASP calculations—full frequency range

Full size image

There are significant changes in the low frequency spectra, that can be more clearly seen in Fig. 11. All VASP results show two low intensity modes at very low frequency with three more intense absorption peaks above 130 cm− 1. Whilst this pattern is the same for all calculations, the actual positions of the peaks vary for the differing methods. The VASP/PAW/GD2 results seem to show the highest shift in frequency from the non-dispersion corrected results with up to 40 cm− 1 shift to higher frequencies in absorption. The GD3, GD3-BJ and TS dispersion correction methods predict absorption spectra in the low frequency regime which are very similar.

Fig. 11
figure 11

IR spectra from VASP calculations—low frequency range

Full size image

The SI gives the full, high, intermediate and low frequency range calculated spectra for all of the methods used. In many respects the observations drawn from the VASP example shown above can be seen in the other methods. There tends to be a small shift to higher frequencies when dispersion corrections are included in the intermediate and high frequency ranges. The intensities are not affected. However in the low frequency range, although the qualitative pattern of absorption is similar, there are significant shifts in the frequency of absorption owing to the inclusion of a dispersion correction. The shift of absorption to higher frequency on the inclusion of dispersion is consistent with the decreased volume of the unit cell, relative to the non-corrected volume. In the cases of CASTEP/NCP19/TS, Crystal/TZVP/GD2 and Crystal/DEF2/GD2 optimization of the S6 parameter resulted in a smaller unit cell and at least in the low frequency regime a shift to higher frequency (see SI).

4.6 Comparison of Calculated Spectra

The calculated spectra were compared with each other by calculating the normalized cross-correlation coefficient between each pair of spectra. This calculation also provides a ‘lag’ or frequency shift which maximizes the cross-correlation for each pair of spectra.

4.6.1 Full Frequency Range Comparison

Figure 12 shows the calculated cross-correlation coefficients for the complete frequency range. The cross-correlation matrix is symmetric and the results presented using a gap-map, where the methods have been clustered and reordered according to their similarity. The clustering is made clear by the dendrogram at the top of the heat-map. The correlation coefficients have been calculated after a lag or shift between each spectrum has been determined which maximizes the correlation coefficient. The heat-map shows the values of the cross-correlation coefficient as a colour map. Yellow is used to describe the highest cross-correlation coefficient (1.0) and blue the lowest (0.5). VASP/PAW, VASP/PAW/GD3 and VASP/PAW/GD3-BJ are shown to be very similar. The none dispersion corrected methods tend to cluster together. Although the CASTEP/NCP19/GD3, CASTEP/NCP19/GD3-BJ and the Abinit GD2 methods also cluster in this region. The GD2 methods form a similar cluster, apart from the Abinit GD2 methods. For CASTEP/NCP19/TS-v the optimization of the S6 coefficient seems to give results which are very similar to the unoptimized result.

Fig. 12
figure 12

Cross-correlation heat-map of full frequency spectra after clustering

Full size image

Figure 13 shows a gap-map created by using the lag frequency to calculate the similarity of each method. The lag frequencies in this plot have been calculated using the full frequency range of the spectra and they vary between − 50 (blue) and + 50 cm− 1 (yellow). Surprisingly this method of clustering shows that the frequency shift which maximizes the correlation between spectra is strongly related to the program used to perform the calculation. It does not seem to be related to the dispersion correction used. All the VASP calculations are clustered together. Crystal/DEF2 and QE/SSSP calculations are clustered together, as are Abinit and CASTEP. Finally all of the Crystal/TZVP calculations are clustered together, probably reflecting the observed trend for Crystal/TZVP calculations to predict lower frequency absorption in the SO4 bending region of the spectrum.

Fig. 13
figure 13

Frequency lag heat-map of full frequency spectra after clustering

Full size image

4.6.2 Low Frequency Range Comparison

The absorption spectrum is dominated by the high frequency region, so it is interesting to see if calculating the cross-correlation function of the low frequency absorption gives similar results. The gap-maps of the correlation function and the lag are provided in the SI. The correlation coefficient shows a wider range of values, ranging from 0.4 to 1.0. Methods with the same dispersion correction seem to be clustered together, as are the non-corrected results. The TS and TS-v results are clustered together in the centre of the table. The least similar groups are the GD2 and non-dispersion corrected results. The GD2 methods (apart from Abinit) cluster together towards the bottom of the Figure. Inspection of the lag heat-map obtained from the low frequency spectrum shows no obvious pattern.

4.7 Comparison of Effective Medium Theories

The results presented so far calculate the effective permittivity of the composite material using a Maxwell-Garnett homogenization formula [9]. Maxwell-Garnett is commonly used in a wide variety of circumstances, but it is not symmetrical with respect to the two components of the composite material. As a result the Bruggeman method [9] is often preferred when the volume fractions of the two components are similar.

For comparison purposes an Averaged Permittivity (AP) effective medium theory [8] with a low volume fraction of sodium peroxodisulfate is used to indicate the position and intensity of absorption from transverse optical phonons with no interaction with the field within the crystal.

Finally, there are occasions where scattering from the particles is important and to understand this calculations have been performed using a Mie methodology, which is relevant for low concentrations of spherical particles embedded in an non-absorbing medium [11].

The calculations presented here use the same parameters as above with only the effective medium method varied. For the purposes of comparing the MG, AP and Bruggeman effective medium methods, Figs. 14, 15 and 16 show the calculated molar absorption spectra from VASP/PAW/GD3-BJ calculations.

Fig. 14
figure 14

Comparison of effective medium theories—VASP/PAW/GD3-BJ high frequency range

Full size image
Fig. 15
figure 15

Comparison of effective medium theories—VASP/PAW/GD3-BJ intermediate frequency range

Full size image
Fig. 16
figure 16

Comparison of effective medium theories—VASP/PAW/GD3-BJ low frequency range

Full size image

Figure 14 shows that in the high frequency range the MG method shifts the absorption to higher frequencies compared with the TO frequencies (shown by AP results), whilst the Bruggeman method produces much broader absorption peaks, with the peak maxima positioned between the AP and MG methods. The intermediate frequency range (Fig. 15) is similar, though the single peak at 550 cm− 1 is independent of the method of calculation. Just above 500 cm− 1 the twin AP peaks are shifted to a slightly higher absorption maximum by the Bruggeman method, whilst the MG method show two peaks both at higher frequency than the AP calculation.

The low frequency spectrum shows the onset of absorption at between 80 and 90 cm− 1. All methods show six absorption peaks. Similarly to the higher frequency ranges, the MG method results in a shift of the absorption maxima to higher frequencies relative to the TO frequencies (shifts of up to 30 cm− 1 are seen), whilst the Bruggeman method tends to show similar, but less marked trends, and much broader absorption peaks. As an example, very similar trends can be seen for the Crystal/DEF2/GD3-BJ calculations reported in the SI.

4.8 Spectra from Mie Scattering Calculations

When the wavelength of light is similar to or smaller than the particles being studied, scattering of light by the particles has to be considered. For spherical particles this can be described well using Mie scattering theory, as long as no multiple scattering events take place. In other words the particles must be very dilute. The spectra shown in Figs. 17 and 18 were obtained from the VASP/PAW/GD3-BJ phonon calculations using the PDielec package. A 10% volume fraction of spheres in PTFE was used with a Lorentzian line width of 5 cm− 1. The smallest particle size (0.1 μm) results for Mie scattering coincide with those of the Maxwell-Garnett method. As the particle size increases the higher frequencies are more affected and the absorption broadens with very little change seen at lower frequencies. Figures showing the full and low frequency ranges can be seen in the SI.

Fig. 17
figure 17

Mie scattering calculations—high frequency range

Full size image
Fig. 18
figure 18

Mie scattering calculations—intermediate frequency range

Full size image

This method can also be used to understand the effect of scattering by air voids or bubbles that are unavoidable in pelletized samples and have shown to contribute to the background spectral response at low frequencies [58, 59].

4.9 Comparison with Experiment

In this section we compare the calculated spectra with experimentally measured IR and THz spectra of Na2(SO4)2. In order to improve the correlation between calculation and experiment and identify any systematic error in the DFT calculations we have also explored re-scaling of the calculated spectra. Such re-scaling is common in molecular calculations where the systematic errors in the calculated frequency of a particular method are corrected by a scale factor [60]. The frequency re-scaling can be expressed as;

$$ f_{\text{new}} = lag + scale f_{\text{calc}} $$
(4)

4.9.1 ATR Spectra

Figure 19 shows a comparison of the experimental ATR spectra with that calculated by VASP/PAW/GD3-BJ using Maxwell-Garnett and Bruggeman effective medium theories with an 80% volume fraction of Na2(SO4)2 in air. For all calculated ATR spectra in this section the effective medium is assumed to be on a slab of diamond with a refractive index of 2.4 with the angle of incidence of the incoming radiation was 45° and the radiation assumed to have equal S and P polarization.

Fig. 19
figure 19

Comparison of experimental ATR signal with Bruggeman and Maxwell-Garnett calculated simulations with scaled frequencies

Full size image

The lag and scale factors used in Fig. 19 are respectively 0.0 cm− 1 and 1.04 for Maxwell-Garnett and 8.6 cm− 1 and 1.04 for Bruggeman.

Both effective medium theories show good agreement with experiment, when re-scaling of the frequency scale is employed. Examination of Table 8 shows that if frequency scaling is used in the comparison of calculated and experimental spectra, all of the methods show a cross-correlation over 0.81 and there is little to choose between the methods. The calculated ATR spectra for all calculations, along with both effective medium approximations for the Crystal/DEF2 calculation can be seen in the SI.

Table 8 scale factors, lag shifts and cross-correlation coefficients between calculated and experimental ATR spectra
Full size table

To compare all the calculated spectra against the experimental spectra a normalized cross-correlation coefficient (as previously discussed in Section 4.6) was calculated between the experimental spectrum in the range 450 to 1400 cm− 1. The calculations were performed with a Maxwell-Garnett effective medium representation of 80% volume fraction of spherical particles Na2(SO4)2 in air. The Lorentzian widths of each transition were chosen so that the peak height of the calculated spectrum agreed with that of the experimental spectrum. The reported cross-correlation coefficients in Table 8 are the maximum coefficients at a constant frequency shift. There are therefore two parameters which are optimized to improve the fit with experiment, a frequency lag and a frequency scale factor. The first three columns in Table 8 show the results for the case that no frequency scaling is employed. The last three columns show the results after optimizing the frequency scaling factor to improve the cross-correlation coefficient.

With no re-scaling of the frequencies the optimum cross-correlation coefficients are found by using a lag shift for the calculated spectra by between 20 and 36 cm− 1 to higher frequency. The combined use of re-scaling and shifting the frequencies results in almost all methods having a cross-correlation with experiment of over 0.8. Only Crystal/TZVP and Crystal/TZVP/GD2 calculations have a cross-correlation below 0.8 and in addition they require a lag shift of over 60 cm− 1. The VASP calculations without re-scaling the frequencies have poor cross-correlations with experiment (below 0.8). However, with re-scaling the cross-correlation is as good as any of the others and the required lag shift in frequency to achieve the best cross-correlation coefficient is small. This behaviour is also shown by CASTEP/NCP19/TS and CASTEP/NCP19/TS-v, where after re-scaling the frequency lag shift required to get the optimum cross-correlation is relatively small. The SI includes a similar comparison for the Bruggeman method. Without re-scaling the frequencies, on average the Bruggeman method requires an additional 9.1 cm− 1 to find the maximum cross-correlation coefficient and the average cross-correlation coefficient is higher by 0.014. However with re-scaling the average cross-correlation coefficient increases by 0.017.

4.10 Transmission Infrared

The experimental transmission infrared spectrum is shown in Fig. 20 and compared with Bruggeman and Maxwell-Garnett effective medium theory calculations based on VASP/PAW/GD3-BJ phonon calculations. The experimental absorption has been re-scaled to show similar peak heights to those calculated and the calculated frequencies re-scaled to improve the position of the calculated peaks. Additional experimental repeat measurements of differing sample concentrations are shown in the SI.

Fig. 20
figure 20

Experimental infrared spectrum for 2.67% mass fraction compared with Bruggeman and Maxwell-Garnett simulations with scaled frequencies

Full size image

In Fig. 20 the values of lag and scale are 4.02 cm− 1 and 1.04 respectively which are similar to those needed for the ATR calculated spectra above with the differences likely owing to the wider and asymmetric peak shapes seen the the experimental spectra. This result indicates a systematic underestimated of the calculated absorption frequencies. Both Bruggeman and the Maxwell-Garnet effective medium theories predict very similar absorption in this region and the same frequency scaling has been applied to both methods.

4.10.1 Terahertz Spectra

Figure 21 shows a comparison of the experimental room temperature terahertz spectrum and the calculated Maxwell-Garnett and Bruggeman effective medium theories using VASP/PAW/GD3-BJ phonon calculations. The experimental spectrum shows a strong background signal which is assumed to arise from scattering of air bubbles trapped in the PTFE supporting matrix [58, 59]. Both simulations reported in the figure account for this scattering through consideration of Mie scattering off a 15% volume fraction of 50 μm air bubbles. The calculated absolute absorption agrees well with experimental measurements. The peak positions and shapes predicted by the Bruggeman effective medium theory are in excellent agreement with experiment. The Maxwell-Garnett method predicts strong absorption at too high a frequency compared with experiment. This shows the choice of effective medium approximation is often crucial at low frequencies to aid in spectral interpretation. Similar results can be seen for a number of the other calculations (not shown) although the best correlation with experiment at these low frequencies is with the VASP/PAW/GD3 phonon calculations using the Bruggeman effective medium approximation.

Fig. 21
figure 21

Comparison of experimental room temperature terahertz measurements with Bruggeman and Maxwell-Garnett effective medium theory simulations

Full size image

5 Conclusions

In this paper we have compared a number of infrared and terahertz spectra of the powdered crystal Na2(SO4)2 to the calculated spectra from a range of DFT programs with various combinations of pseudo-potentials and van der Waals’ dispersive corrections. The inclusion of a van der Waals’ dispersion correction has a significant effect on the calculated absorption spectrum, and are crucial for good correlation at low frequencies, where there is a systematic shift of absorption to higher frequencies; shifts of over 40 cm− 1 were seen. The default values of the dispersion correction parameters S6 and SR have been determined for a wide range of molecules and optimizing these parameters to improve the predictions for a single molecule can lead to poor results, especially if only a single parameter such as the volume is chosen for improvement. Determining an optimum parameter does have an impact on the predicted spectrum, but generally speaking it is smaller than other factors in the calculation.

Low frequencies were particularly influenced by aliasing issues associated with the grids used to store the wavefunction and charge. But in all cases projection of the crystal translation from the dynamical matrix provided sensible results.

For the plane-wave based calculations convergence of the calculation was relatively straightforward to achieve. Although it is important to confirm that properties such as the phonon frequencies have converged as well as the lattice energy. For Crystal it was more difficult to be sure that the atom centred basis set was adequate. The TZVP basis set was not adequate and the DEF2 basis which was slightly larger gave similar results to the plane-wave calculations.

The use of the cross-correlation of the predicted spectra to generate gap-maps of the similarities between the calculation was a useful tool and highlighted how the TZVP basis stood out from the other calculations. It also showed how the use of cell volume to determining S6 in Crystal/DEF2/GD2-v resulted in an absorption spectrum which was different to other calculations.

In the calculations of the absorption spectrum, the use of effective medium methods to calculate the changes in absorption frequency and intensity owing to the interaction of the electromagnetic radiation field with the internal field generated by the vibrating crystal, was important over the whole frequency range. The Maxwell-Garnett method predicted the largest changes with shifts to higher frequencies of up to 40 cm− 1. The Bruggeman method tends to show broader absorption at frequencies intermediate between the TO frequencies and the absorption maxima predicted by Maxwell-Garnett. For small particle sizes relative to the wavelength of the radiation, the Mie method for incorporating scattering effects agrees well with the Maxwell-Garnett effective medium theory. For particles smaller the 1 μm the Maxwell-Garnett and Mie methods agree well up to 600 cm− 1, above this frequency the transitions get broader and show additional scattering artifacts in the calculated absorption. The incorporation of Mie scattering from air bubbles trapped in the support matrix greatly improves the agreement of the calculated THz spectrum with the experimental spectrum.

All the post analysis methods described in this paper including the effective medium approximations, Mie scattering, cross-correlation with experiment and the optimization of lag shift and scale factors are available in the latest release of PDielec [61].