DFT Calculations of Temperature-Dependent NQR Parameters in {\alpha}-paradichlorobenzene and {\beta}-HMX

TL;DR

This paper presents a first-principles DFT-based method to predict temperature-dependent NQR spectra in molecular crystals without relying on empirical potentials or detailed internal motion models, applicable to any periodic solid with known structure.

Contribution

The method combines DFT, molecular dynamics, and dynamic corrections to accurately predict T-dependent NQR parameters ab initio, independent of specific intermolecular potential assumptions.

Findings

Successfully predicts temperature-dependent NQR frequencies for 35Cl and 14N in two crystals.

Requires no prior knowledge of intermolecular potentials or internal motions.

Applicable to any periodic solid with known structural data.

Abstract

A method for first principles predictions of observed temperature-dependent NQR spectra is presented using density functional theory (DFT) and the isobaric T-dependent NQR frequencies of 35Cl and 14N nuclei are computed for the two molecular crystals (1) alpha-paradichlorobenzene, and (2) the nitroamine high explosive beta-HMX (beta-octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine) over a range of 200 K and up to room temperature. Notably, the method requires no supposition of a form for the intermolecular potentials or wave functions, requires no particular insight as to the nature of the internal motions or of the chemical bonds present, and does not depend on the crystal structure, making the method amenable to any periodic solid for which experimentally determined structural data are available. For each substance, unit cells of various volume are prepared using experimentally determined atomic positions and cell parameters. In each of the prepared volume-corrected cells, a molecular dynamics (MD) simulation generates a set of perturbed atomic positions, the collection of which is intended to represent the system at a given T,V. For each configuration of atoms generated along the MD trajectory, the electric field gradient (EFG) tensors are computed at the site of each quadrupolar nucleus. The rotational displacements of the moving EFG principal axes from their equilibrium directions are used to apply a dynamic correction to the DFT-computed static-lattice NQR frequencies, resulting in first-principles DFT predictions of T-dependent NQR spectra at constant pressure. Because the V-dependence from thermal expansion and the T-dependence due to internal motions are simultaneously considered, the NQR's notoriously model-dependent temperature coefficients are computed entirely ab-initio.

Figures

Tables

Table 1. Table 1: The experimental cell parameters of α 𝛼 \alpha -pCl ϕ 2 subscript italic-ϕ 2 {}_{2}\phi Wheeler and Colson ( 1976 ) Estop  et al. ( 1997 ) and β 𝛽 \beta -HMX Deschamps  et al. ( 2011 ) used in this study. The the variation of the lattice constants can, in general materials, either increase or decrease C q , 0 subscript 𝐶 𝑞 0 C_{q,0} . In molecular crystals, increasing volume typically increases C q , 0 subscript 𝐶 𝑞 0 C_{q,0} , as intermolecular separation leads to greater anisotropy of the nuclear electronic environment. Our results in table 2 agree with this trend, making volume and temperature have competing effects on the observed NQR frequency, though the volume dependence of C q , 0 subscript 𝐶 𝑞 0 C_{q,0} was found to be nearly an order of magnitude stronger in α 𝛼 \alpha -pCl ϕ 2 subscript italic-ϕ 2 {}_{2}\phi than in β 𝛽 \beta -HMX as shown in figure 5 .

	T (K)	a $\left(\text{Å}\right)$	b $\left(\text{Å}\right)$	c $\left(\text{Å}\right)$	$\beta$	Volume $\left({\text{Å}}^3\right)$
	293	14.754	5.840	4.025	112.52∘	320.3
	273	14.747	5.830	4.014	112.33∘	319.2
$\alpha$-pCl${}_2{}^{}\varphi$	250	14.730	5.812	3.997	112.24∘	316.7
	225	14.720	5.801	3.982	112.17∘	314.8
	200	14.705	5.787	3.967	112.00∘	311.4
	100	14.664	5.740	3.925	111.77∘	306.8
	303	6.5255	11.0369	7.3640	102.670∘	517.45
	293	6.5245	11.0240	7.3619	102.642∘	516.68
	273	6.5289	10.9875	7.3453	102.616∘	514.29
	248	6.5254	10.9702	7.3503	102.582∘	513.53
$\beta$-HMX	223	6.5334	10.9419	7.3421	102.491∘	512.45
	198	6.5206	10.9123	7.3395	102.467∘	509.93
	173	6.5273	10.8834	7.3286	102.366∘	508.54
	123	6.533	10.8400	7.3207	102.271∘	505.97

Table 2. Table 2: Comparison of the computed T-dependent NQR parameters from this work with experimental values in α 𝛼 \alpha -pCl ϕ 2 subscript italic-ϕ 2 {}_{2}\phi Moross and Story ( 1966 ) and β 𝛽 \beta -HMX. Landers  et al. ( 1981 ) The parameter C q , 0 subscript 𝐶 𝑞 0 C_{q,0} is derived as an average of the C q , 0 i superscript subscript 𝐶 𝑞 0 𝑖 C_{q,0}^{i} computed with the i t ​ h superscript 𝑖 𝑡 ℎ i^{th} output of the molecular dynamics simulation. κ T subscript 𝜅 𝑇 \kappa_{T} is given by equation ( 12 ). Note that mean values of C q , 0 = 6.11 subscript 𝐶 𝑞 0 6.11 C_{q,0}=6.11 MHz have been imputed for HMX’s equatorial sites at 273K and 293K because the computed averages showed signs of systematic error. The absolute error column contains the absolute value of the difference between the calculated and experimental NQR frequency ν ​ ( T ) 𝜈 𝑇 \nu(T) .

		computed NQR parameters					experimental	error
structure	T (K)	$C_{q,0}$ (MHz)	${\eta }_0$	${\nu }_0$ (MHz)	${\kappa }_T$	$\nu \left(T\right)$ (MHz)	$\nu \left(T\right)$ (MHz)	$\Delta \nu$ (MHz)
	293	69.3100	0.09749	34.7099	0.9729	33.7696	34.2869	0.5486
	273	69.3241	0.09714	34.7165	0.9728	33.7725	34.3505	0.5424
$\alpha$-pCl${}_2{}^{}\varphi$	250	69.2976	0.09687	34.7029	0.9788	33.9679	34.4190	0.5771
	225	69.2819	0.09650	34.6947	0.9774	33.9109	34.4880	0.4511
	200	69.2680	0.09635	34.6875	0.9804	34.0088	34.5511	0.5780
	100	69.2026	0.09566	34.6540	0.9868	34.1966	34.7452	0.5173
	303	6.2343	0.5312	5.5036	0.9852	5.4219	5.3008	0.1109
	293	6.2147	0.5284	5.4821	0.9891	5.4226	5.3026	0.1105
	273	6.2070	0.5295	5.4769	0.9900	5.4219	5.3061	0.1221
$\beta$-HMX	248	6.2117	0.5312	5.4837	0.9895	5.4262	5.3105	0.1195
(axial sites)	223	6.2168	0.5313	5.4884	0.9902	5.4344	5.3149	0.1156
	198	6.2133	0.5302	5.4835	0.9923	5.4413	5.3191	0.1157
	173	6.2077	0.5318	5.4811	0.9914	5.4337	5.3233	0.1200
	303	6.1314	0.4850	5.3420	0.9847	5.2601	5.0643	0.1984
	293	6.1100	0.4833	5.2642	0.9871	5.2521	5.0657	0.1964
	273	6.1100	0.4833	5.2593	0.9883	5.2587	5.0685	0.1945
$\beta$-HMX	248	6.1152	0.4813	5.3222	0.9901	5.2693	5.0718	0.2026
(equatorial sites)	223	6.1225	0.4832	5.3314	0.9899	5.2779	5.0751	0.1974
	198	6.1052	0.4830	5.3160	0.9919	5.2729	5.0784	0.1903
	173	6.1096	0.4838	5.3211	0.9919	5.2778	5.0815	0.1864
	123	6.1002	0.4838	5.3130	0.9948	5.2854	5.0869	0.1957

Table 3. Table 3: Fitting parameters of calculated and experimental NQR frequencies for functional form ν ​ ( T ) = a 0 + a 1 ​ T + a 2 ​ T 2 + a − 1 / T 𝜈 𝑇 subscript 𝑎 0 subscript 𝑎 1 𝑇 subscript 𝑎 2 superscript 𝑇 2 subscript 𝑎 1 𝑇 \nu\left(T\right)=a_{0}+a_{1}T+a_{2}T^{2}+a_{-1}/T .

structure	dataset	fitting parameters
		$a_0$ (MHz)	$a_1$ $\left(\frac{\text{MHz}\cdot {\text{K}}^{-1}}{{10}^{-4}}\right)$	$a_2$ $\left(\frac{\text{MHz}\cdot {\text{K}}^{-2}}{{10}^{-6}}\right)$	$a_{-1}\left(\text{MHz}\cdot \text{K}\right)$
	experimental	34.8437	-5.2424	-4.6970	—
$\alpha$-pCl2-$\varphi$	calculated	34.2932	-5.5595	-4.2895	—
	error (no units)	1.56%	25.8%	12.9%	—
	experimental	5.3761	-2.2642	—	-2.3392
$\beta$-HMX (axial)	calculated	5.4885	-1.9402	—	-2.8598
	error (no units)	2.09%	14.3%	—	22.3%
	experimental	5.1372	-2.0976	—	-3.3160
$\beta$-HMX (equatorial)	calculated	5.3426	-2.5139	—	-3.3235
	error (no units)	4.00%	19.84%	—	0.226%

Equations

$$(\nabla \mathbf{E}{)}^{ij}=\frac{{\partial }^2\varphi }{\partial x_i\partial x_j}={\varphi }_{ij}$$

$$C_q=\frac{eQ{\varphi }_{ZZ}}{h}$$

$$\eta =\frac{{\varphi }_{YY}-{\varphi }_{XX}}{{\varphi }_{ZZ}}$$

$$\nu =\frac{1}{2}{C}_q\sqrt{1+\frac{{\eta }^2}{3}}$$

$${\nu }_{\pm }=\frac{3}{4}{C}_q\left(1\pm \frac{\eta }{3}\right)$$

$${\nu }_{\delta }={\nu }_{+}-{\nu }_{-}=\frac{\eta C_q}{2}$$

$${\nu }_q\left(T\right)=a-\frac{b}{\exp \left(c/T\right)-1}$$

$${\nu }_q\left(T\right)={\nu }_0\left(1+b^{\prime }T+\frac{c^{\prime }}{T}\right)=a+bT+\frac{c}{T}$$

$${\varphi }_{X^{\prime }{X}^{\prime }}$$

$${\varphi }_{Y^{\prime }{Y}^{\prime }}$$

$${\varphi }_{Z^{\prime }{Z}^{\prime }}$$

$${\varphi }_{X^{\prime }{Y}^{\prime }}$$

$${\varphi }_{Y^{\prime }{Z}^{\prime }}$$

$${\varphi }_{Z^{\prime }{X}^{\prime }}$$

$$⟨{\theta }_i⟩=0$$

$$C_q\left(T\right)={\kappa }_T{C}_{q,0}$$

$$\kappa_{T}=\Big{(}1-\frac{3}{2}\left(\langle\theta_{X}^{2}\rangle+\langle\theta_{Y}^{2}\rangle\right)-\frac{1}{2}\eta\left(\langle\theta_{X}^{2}\rangle-\langle\theta_{Y}^{2}\rangle\right)\Big{)}$$

$$\nu \left(T\right)={\kappa }_T{\nu }_0$$

$$\nu \left(T\right)$$

$$=a_0+a_{1}T+a_{-1}/T$$

$$\nu \left(T\right)$$

$$=a_1^{\prime \prime }T\left(1+a_1^{\prime }T+a_{-1}^{\prime }/T\right)$$

$$=a_0+a_{1}T+a_2{T}^2$$

Full text

nuclear quadrupole resonance, density functional theory, gauge including projector augmented wave

I Introduction

Nuclear quadrupole resonance (NQR) spectroscopy is a radio-frequency (RF) technique for solid materials in which a quadrupolar nucleus is irradiated with RF magnetic fields and the resulting transitions between the nuclear spin energy levels are directly observed. Pound (1950) NQR arises from the interaction between a nuclear electric quadrupole moment and the local electric field gradient (EFG) that is natural to the system. Because the EFG, a second-rank tensor, depends on the particular charge distribution in the space surrounding the nucleus, NQR is a sensitive probe of the electronic structure of systems. As such, NQR experiments have long been used to carry out basic investigations on the structure and dynamics of a wide range of solid materials. Dodgen and Ragle (1956) Kind and Roos (1976) Karpowicz and Brill (1983) Karpowicz and Brill (1984)

Calculations of NQR parameters can assist in the interpretation of NQR spectra and help overcome practical difficulties of carrying out NQR measurements. The inability to predict NQR unknown frequencies presents the greatest challenge for NQR experimentation in new materials. This study presents a method to overcome this challenge by allowing direct calculations of room temperature NQR frequencies that are sufficiently accurate to narrow the search space for NQR exploration in new materials.

I.1 NQR Applications and Limitations

Because NQR frequencies in a material result from the substance’s unique electronic structure, detection of the material-specific NQR spectra can unambiguously indicate the presence of a particular substance Das and Hahn (1958) and can be used to identify or differentiate materials. Exploitation of ths ”chemical fingerprint” property of NQR spectra has long been considered in the development of technology of considerable practical use. NQR is known to be capable of identifying counterfeit pharmaceuticals Dass and Gajjar (2011). Additionally, NQR has a long history of consideration for remote and non-invasive detection of materials such as landmines, bombs, Buess et al. (1993) Smith (1995) Mozjoukhine (2000) Gregorovič and Apih (2009) and narcotics Garroway et al. (1994) Magnuson et al. (2001), as is evidenced by the patent record of the past three decades. The United States government has extensively characterized NQR spectra of energetic materials like TNT and HMX. Marino et al. (1982) Yet, despite the obvious merits of a point-and-shoot bomb detector, and the substantial body of NQR research conducted to this end, no such device has been shown to work in practice. Nevertheless, development of NQR continues to this day in the pursuit of miniaturization and automation of detector operation, including ”smart” devices running system-on-chip operating systems and other technologies Ariando et al. (2019).

While there is continuous development of NQR technology for its potential commercial use, NQR is much less common than NMR Vij (2007) studies in basic solid state research due to practical difficulties of the technique. For the inductive pickups used in NQR and NMR detection, signal-to-noise ratio (SNR) increases with frequency due to Faraday’s law of induction. Yet the naturally arising NQR frequencies are often quite low: most NQR frequencies fall below 100MHz but can be as low as several kHz Lucken (1969) in the case of ${}^{14}N$ due to its small quadrupole moment. Consequently, NQR studies can have prohibitively poor (SNR) Das and Hahn (1958). NMR studies often have much higher SNR because the Larmor frequency $\omega_{L}=\gamma B_{ext}$, and hence the SNR, can be made very high by increasing the applied magnetic field $B_{ext}$, where $\gamma$ is the gyromagnetic ratio of nucleus.

Additionally, NQR is often shelved in favor of NMR because typical NQR detectors are easily reused as NMR detectors. Because they are determined by the system’s natural EFG, which is a consequence of the system’s particular electronic structure, NQR frequencies are material specific; observed NQR frequencies cannot be chosen or appreciably modified by an experimenter. As such, the operating frequency of an NQR detector is determined by the sample itself. The low SNR in NQR studies necessitates the receiver section of an NQR detector to have very high sensitivity in a narrow passband to accommodate typical NQR line-widths of a few kHz. Due to the impedance matching needs of the sample coil pickup circuitry, the consequences of lossy substrates, and other challenges in the construction of necessary low noise amplifiers and probe circuits NQR detectors that are capable of broadband operation is a challenging engineering task. As a consequence, NQR experiments are often carried out using home-built ”one-off” spectrometers, useful for only a sample under consideration, or at best, some other similar materials. Development of off-the-shelf, reusable NQR spectrometers remains an active area of research for commercial use, however. Dass and Gajjar (2011)

The principal reason NQR methods are not used as extensively as NMR for investigating the properties of materials is the inability to predict the frequencies accurately. Even approximate theoretical prediction of NQR frequencies would make NQR easier to carry out. This is because for new materials, in which NQR has never been performed, locating unknown NQR requires a search over a range of plausible NQR frequencies as they are not known a-priori to the experimenter. Low SNR and ignorance of the NQR frequencies make NQR impractical or impossible if the search space is not sufficiently constrained, as one may not observe the NQR frequencies at all. Accurate calculations of unknown NQR frequencies by computer simulation would mitigate the problem of locating unknown NQR lines, reducing the barrier to entry for NQR use to probe interesting physics in new materials.

Recent developments in DFT codes allow calculation of NQR frequencies that are valid for $T=0$ K. Bonhomme et al. (2010) However NQR experiments (or any experiment, for that matter) are more easily performed at room temperature, making reliable prediction of room temperature NQR frequencies highly desirable. A method for making such predictions with a straight-forward, broadly applicable procedure for use in arbitrary materials is presented here.

I.2 GIPAW NQR Predictions and Limitations

Introduced in 2001 to compute NMR chemical shifts, the gauge-including projector augmented-wave (GIPAW) method Pickard and Mauri (2001) allows calculations of electric field gradients and corresponding static-lattice NQR frequencies in the plane wave pseudopotential formalism of density-functional theory (DFT) Yates and Pickard (2008) Bonhomme et al. (2010) Charpentier (2011) Bonhomme et al. (2012). However, the strong temperature dependence of NQR frequencies, which arises primarily due to internal motions, can not be ascertained using the vanilla DFT routines that are currently available because the GIPAW method, like many routines based on DFT, can only be applied to a static-lattice undergoing no internal motions. Without consideration of lattice dynamics, the results of GIPAW calculated NQR frequencies are only valid for $T=0$ K. Yet, all the practical commercial uses of NQR detection in materials detection and characterization would seem to require room temperature operation; after all, airports tend to be kept at room temperature and pharmaceutical use cases would be more realistic without the need for cryogenics. Even for basic science needs, NQR as a probe of interesting transitions across a temperature range involves slowly cooling down the sample from higher temperatures. Because NQR frequencies can vary hundreds of kHz or more across a temperature range of 300 K, static lattice calculations do not sufficiently narrow the search space for finding unknown NQR frequencies.

GIPAW NQR calculations that are performed in the static lattice necessarily neglect lattice dynamics, so should not be expected to match room temperature NQR frequencies. Furthermore, unless the atomic positions and lattice vectors used in the GIPAW calculation were taken at the same temperature as the NQR experiment, the often significant volume dependence of the NQR frequencies is also neglected by the calculations.

Presently, no GIPAW study known to this author accounts for the strong temperature dependence of NQR frequencies. Previous work on the development of GIPAW calculations of NQR frequencies is often performed in the static-lattice in a variety of structures. Socha et al. (2017) Milinković and Bilalbegović (2012) Perić et al. (2014) Leroy et al. (2019) The calculations not only neglect the effect of motions, but furthermore, if the simulation’s cell volume does not match that of the NQR experiment’s, then the calculations are performed on a system that is perhaps very different from the experimental one. The relevance of these calculations to room temperature NQR data is unclear, though the comparison is often made. Thermal expansion can have a number of effects on NQR frequencies. If the structure is not optimized with the correct lattice constants, a DFT simulation would not faithfully represent the “experimental static-lattice NQR frequency“ - that is, the NQR one would measure in a hypothetical frozen lattice that has the cell volume of the room temperature sample. Not even the fitting parameter meant by “static-lattice NQR frequency“, which is never directly observed, is computed as advertised if care is not taken to prepare the cell to faithfully represent the structure studied in the experiment.

I.3 Inclusion of Motional Effects

Consideration of lattice dynamics in conjunction with DFT calculations is indeed a fundamental problem. The effect of internal motions on DFT calculated NMR shielding tensors has been investigated using vibrational averaging in both solid materials Lee et al. (2007) Rossano et al. (2005) as well as gases and solutions Tang and Case (2007) Bagno et al. (2007) Waller et al. (2008) Röhrig and Sebastiani (2008) Kongsted et al. (2007) Straka et al. (2008) Harding et al. (2008). Molecular dynamics (MD) simulations have been used in conjunction with GIPAW to include the influence of internal motions in DFT calculated NMR chemical shifts in organic solid materials Dumez and Pickard (2009), but no method has been proposed to compute temperature dependent NQR parameters in materials from first principles, but no method has been proposed to account for the effect of internal motions on observed NQR using first principles DFT calculations. In this study, the plane-wave DFT code QUANTUM ESPRESSO Giannozzi et al. (2009) Giannozzi et al. (2017) is used to determine the isobaric $T$-dependence of 35Cl and 14N NQR frequencies in $\alpha$-pCl${}_{2}\phi$ and $\beta$-HMX, respectively.

By accounting for the influence of internal motions and the volume dependence of the electronic structure simultaneously, the T-dependent NQR spectra and corresponding BBK model temperature coefficients are predicted at atmospheric pressure in both materials using the same procedure. Notably, this method does not require any consideration of the vibrational spectra, nor of the nature of the structure and chemical bonding in the system. Taking as inputs only the relevant crystallographic structural information of atomic positions and the effect of thermal expansion on lattice vectors, the method is both structure-agnostic and model-independent making its application to arbitrary systems straight-forward.

II Theory

II.1 NQR Transition Frequencies in the Static Lattice

Nuclei with a non-spherical nuclear charge density $\rho(\mathbf{x})$ can have an electric quadrupole moment

${Q=\frac{1}{e}\int(3z^{2}-r^{2})\rho(\mathbf{x})d^{3}x}$ that interacts with the external electric potential $\phi$ due to its electronic environment via the electric field gradient (EFG)

[TABLE]

Similar to the manner of NMR, in which the interaction of the nuclear magnetic dipole moment with an external magnetic field can be probed by irradiation of the nucleus with RF of a particular frequency, NQR frequencies result from the splitting of the nuclear spin energy levels caused by this interaction. NQR frequencies depend on the quadrupole coupling constant $C_{q}$ and the asymmetry parameter $\eta$ defined as

[TABLE]

[TABLE]

where $\phi_{ii}$ are the eigenvalues of the EFG tensor at the nuclear site. The directions of the eigenvectors define the principle axes of the EFG tensor, the names of which are chosen by the convention $|\phi_{ZZ}|\geq|\phi_{YY}|\geq|\phi_{XX}|$ ensuring $0\leq\eta\leq 1$. NQR transition frequencies are then expressed in terms of only $C_{q}$ and $\eta$, though the form of the expression will depend on the nuclear spin of the particular isotope. For spin $I=\frac{1}{2}$, $Q=0$ and there are no transition frequencies. If $I>\frac{1}{2}$, the NQR frequencies are proportional to $C_{q}$ and a factor in $\eta$ which can have a significant impact on $\nu_{q}$ when $\eta$ approaches unity.

35Cl has nuclear spin $I=\frac{3}{2}$ in which case there is just one NQR frequency

[TABLE]

In the $I=1$ case of 14N, there can be one, two, or three transitions for $I=1$ depending on the magnitude of $\eta$. The frequencies are given by

[TABLE]

[TABLE]

Note that as $\eta$ approaches zero, $\nu_{+}\to\nu_{-}$ and $\nu_{\delta}\to 0$. As $\eta$ approaches unity, $\nu_{\delta}\to\nu_{-}=\frac{1}{2}\nu_{+}$.

II.2 Motional Effects and T-dependence

As a nucleus is displaced due to motions within a crystal, it is subjected to a rapidly changing electronic environment. The nuclear position changes at a high frequency within the spatially dependent EFG. Furthermore, EFG is also time dependent in the constantly evolving system undergoing motions. The EFG seen by the nucleus therefore has a somewhat intractable dependence on space and time. Yet observed NQR frequencies are not time dependent because the motional frequencies in solids are many orders of magnitude higher than the NQR transition frequencies. The effective observed $C_{q}$ in a system with a dynamic lattice is determined by the time averaged EFG to which the nucleus is subjected in a system undergoing motion. The nature of the internal motions change significantly with temperature as the magnitude of the nuclear displacements and number of modes increase with with increasing $T$. The internal motions that lead to the time averaged effective EFG at the nuclear site are the cause of NQR’s strict $T$-dependence. The time averaged EFG can be shown to have a weaker associated $C_{q}$ than that of the equilibrium structure, meaning internal motions act to decrease the observed NQR frequency.

While the effect of motions is normally to reduce the NQR frequencies, the nature of the $V$-dependence due to thermal expansion (an implicit $T$-dependence) is more subtle. If the expansion is highly anisotropic, as is the case in many monoclinic crystals, both $C_{q}$ and $\eta$ are affected by changes in the magnitudes and directions of the lattice vectors. Thermal expansion can also cause reorientation of molecules within the cells of molecular crystals, which may affect the NQR parameters in any fashion, either gradually or abruptly. Moreover, in molecular crystals, intermolecular distances increase during expansion of the cell, while intramolecular atomic distances do not appreciably change. As a result, increasing molecular separation in crystal cells gives rise to EFG’s with larger associated magnitudes of $C_{q}$ as the electronic environment of the nucleus possesses greater anisotropy in a lone molecule than in a crystal cell with neighboring molecules. This means that in molecular crystals, $C_{q}$ and the NQR frequency $\nu_{q}$ at constant $T$ strictly increase with increasing $V$. On the contrary, ionic crystals display the opposite trend because the limiting case as $V\to\infty$ is for the nucleus to be alone in space, making the local EFG vanish. Because of the complexity of NQR $V$-dependence, previous models for the precise nature of the isobaric NQR $T$-dependence vary between structures as the volume correction to observed EFGs can in general have any affect on $\nu_{q}$. In this work, a method that is generally applicable to materials irrespective of their nature is proposed to account for all these effects. By performing a particular sequence of DFT calculations resulting in the calculation of the correct $C_{q}$ and $\eta$, along with a dynamic normalization factor $\kappa_{T}$ which accounts for internal motions, the computed parameters are sensitive to both $T$ and $V$, allowing one to predict laboratory NQR $T$-dependence.

II.3 Pioneering Work

A general observation from early NQR research is that the rate of change of the NQR frequency with respect to temperature $\frac{d\nu}{dT}$ tends to zero at low temperature and becomes a negative constant as $T$ increases to 293 K and beyond. Brown (1990). It was found that the character of the $T$-dependence differed in ionic versus molecular crystals. Furthermore, the values of $C_{q}$ in a molecular crystal were typically around 5% lower than the ones observed in single molecules. Attribution of these effects to underlying physical processes, disentangling their dependence on $T$ and $V$, accounting for otherwise anomalous behavior of $T$-dependent NQR frequencies became the subject of many theoretical works in the following decades which carefully consider the dynamics and the thermodynamic properties of individual systems.

The Bayer model is a theoretical framework for the $T$-dependence of NQR frequencies that was first reported by Horst Bayer in 1951 Bayer (1951). The concepts in the paper underpin many similar models which are currently in use that are essentially extensions to his work accounting for additional effects. Bayer proposed that the internal motions of the system are responsible for the temperature curves observed in NQR experiments by considering the rotations of the EFG principle axes in a system undergoing motion. Paradichlorobenzene was used to illustrate this concept in his original paper and it is shown that internal motions result in a reduced NQR frequency compared with that which would be observed in a static lattice. Bayer arrived at a functional form for NQR frequencies which captures the desired limiting behavior

[TABLE]

where $a$, $b$, and $c$ are parameters to be fit. The model considers a single model of librations for the paradichlorobenze molecule. Bayer’s theory had some success but because it is a constant volume theory, it failed to account for the variation of $\nu_{q}$ across larger temperature ranges if thermal expansion significantly impacts the static lattice EFG.

Highly influential experimental and theoretical work by Kushida, Benedeck, Bloombergen Kushida et al. (1956) reported in 1956 presented KBB model which extends the ideas of Bayer to include arbitrary numbers of vibrational modes and the effect thermal expansion on the parameter $C_{q}$ 111In the KBB paper as well as many works, the parameter $C_{q}$ is referred to instead by $q=\phi_{ZZ}/e$ and the resulting NQR frequency. The KBB paper emphasizes the thermodynamic considerations affecting NQR frequencies by separately considering the dependence of NQR parameters temperature, pressure, and volume. The KBB paper reported very thorough measurements of the pressure dependence of NQR frequencies in three polymorphs of paradichlorobenzene at several temperatures. Their work elucidated the complexity of the general problem of modeling NQR’s isobaric $T$-dependence by disentangling the various competing effects of each $T,V,P$, as well as distinguishing ionic and molecular crystals by including measurements $P$-dependence NQR frequencies at constant T in potassium chlorate (KClO3) and cuprous oxide (Cu2O). The KBB model proposes a functional form

[TABLE]

where $\nu_{0}=a$ refers to the NQR frequency computed by equations (4) and (5) using the values of $C_{q}$ and $\eta$ at $T=0$. KBB acknowledges the volume dependence of $C_{q}$ and $\eta$, which the authors addressed by assuming an equation of state for the solid in order to construct $\nu$-$vs$-$V$ plots. Their work emphasized the relationship of NQR parameters to thermodynamic conditions, characterizing the individual dependence on $P,V$ and $T$ of observed NQR frequencies in several materials.

The effect on $\nu_{q}$ due to the $V$-dependence of the amplitudes of lattice vibrations was carried out Brown (1960) shortly after giving rise to the Brown model in which a quadratic term in $T^{2}$ is added to the KBB fitting function. The seminal works of Bayer, KBB, and Brown have stood the test of time as their early considerations underpin many existing variants of their original functional forms for $\nu_{q}(T)$, which primarily offer corrections for additional effects and special cases. A summary of just some of the many available fitting functions is presented in one of the early attempts to study NQR frequencies using DFT by simulating the vibrational spectra of molecules in order to directly calculate the fitting parameters for the Bayer, KBB, and Brown models, as well some of their many descendants. Latosińska et al. (2002)

II.4 The Dynamic Normalization Factor $\kappa_{T}$

In this work, the many effects considered separately by the Bayer, KBB, and Brown models are captured at once by determination of a parameter $\kappa_{T}$ which accounts for internal motions, and by explicit calculation the $V$-dependence of $C_{q}$ and $\eta$. Calculation of the isobaric $T$-dependence of NQR frequencies amounts to calculation of $\kappa_{T}$, $C_{q}(V)$, and $\eta(V)$. It is worthwhile to review the original assumptions of Bayer which underpin the theoretical basis of this work. Clear descriptions of these essential considerations were given by Hahn Das and Hahn (1958) and Lucken shortly later Lucken (1969).

Consider a nucleus that is displaced from its equilibrium position in a system undergoing motion. The displaced nucleus will experience at its new location an EFG tensor with principle axes $X$, $Y$, and $Z$ having orientations that are rotated with respect to the directions of the equilibrium principle axes $X^{\prime}$, $Y^{\prime}$, $Z^{\prime}$.

The rotated tensor components can be expressed in the space-fixed coordinate frame as Das and Hahn (1958) Lucken (1969)

[TABLE]

where $\theta_{i}$ are intrinsic rotations about the $i^{th}$ axis through which the stationary equilibrium principle axes system could be mapped to the moving axes. It is proposed that the observed coupling constant for a system at non-zero temperature will be given by an average of the $\phi_{i^{\prime}j^{\prime}}$ in equation (II.4) over the time-dependent angular displacements $\theta_{X},\theta_{Y}$ and $\theta_{Z}$ Lucken (1969) the system develops as a result of its various modes internal motion. Presuming the motions are a superposition of simple oscillators, the average angular displacements vanish,

[TABLE]

averages of the squared angular displacements will survive in equation (II.4). Therefore the averaged tensor is diagonal and it can be shown that the $T$-dependent coupling constant is given by the static lattice value of the coupling constant $C_{q,0}$ and a dynamic normalization factor $\kappa_{T}$, that is depends on both $V$ and $T$, which accounts for the effect internal motions

[TABLE]

where $C_{q,0}$ is the volume ($V$) dependent static-lattice value of the coupling constant$C_{q}$ and $\kappa_{T}$ is given by Wang (1955)

[TABLE]

The NQR frequencies are always proportional to $C_{q}$ and therefore are adjusted by the same factor $\kappa_{T}$,

[TABLE]

where $\nu_{0}$ is the NQR frequency calculated using the static lattice formulas (4) or (5) and $C_{q,0}$.

$C_{q,0}$ and $\eta$ will, of course, generally depend on volume in a complex manner due to overall changes in the structure and orientation of the system and it?s constituent nuclei or molecules, was noted in section II.2. The subtlety of volume dependence and its net effect on NQR frequencies is the primary roadblock in to developing a model for isobaric $T$-dependence that applies to general materials and temperature ranges. Equations (11) and (13) imply that the isobaric $T$-dependent NQR frequencies can be computed exactly if $C_{q,0}$, $\eta$, and $\kappa_{T}$, and their dependence on $T,V$ are known.

Because $\langle\theta_{i}^{2}\rangle>0$, and $\kappa_{T}$ is reduced by $\frac{3}{2}\left(\langle\theta_{X}^{2}\rangle+\langle\theta_{Y}^{2}\rangle\right)$, which dominates the term in $\eta$ and $\langle\theta_{X}^{2}\rangle-\langle\theta_{Y}^{2}\rangle$, $\kappa_{T}$ decreases as the magnitude of angular displacements increases. Generally, higher temperature leads to more energetic internal motions with larger values of $\langle\theta_{i}^{2}\rangle$, so $\kappa_{T}$ should decrease with increasing $T$. The observed NQR frequency $\nu\left(T\right)=\kappa_{T}\nu_{0}$ will then also decrease with increasing $T$ if $\nu_{0}$ does not increase with volume faster than $\kappa_{T}$.

A net decrease in NQR frequency with increasing $T$ at atmospheric pressure is observed in most materials, though counter-examples exist Kind and Roos (1976) when either gradual or abrupt structural changes increase $\nu_{0}$ so sharply that the affect of $\kappa_{T}$ is overcome. On the other hand, the volume variation of $\nu_{0}$ may sometimes be ignored, supposing that $\nu_{0}$ is a constant. Whatever the dependence of $\nu_{0}$ on volume in a given material in fact is, $\nu_{0}(V)$ is readily accessible obtain using GIPAW calculations: simply prepare a set of crystal cells with lattice vectors and atomic positions each commensurate with the desired $T,V$, then, compute the static lattice NQR parameters $C_{q,0}$ and $\eta$ in each cell.

III Methods

NQR frequencies depend on $T$ directly through internal motions and also depend on volume (and so $T$) as thermal expansion of the material causes structural changes that modulate all of the NQR parameters in potentially subtle ways. In this study, both effects are accounted for by straight-forward computation of the volume dependence of $C_{q,0}$ and $\eta$, and the dependence of $\kappa_{T}$ on both $T$ and $V$. The NQR frequency at a given $T,V$ is completely determined by $C_{q,0}$, $\eta$, and $\kappa_{T}$ by equations (4), (5), and (13). Accurate calculation of $\kappa_{T}$, $C_{q,0}$, and $\eta$ by a sequence of individual DFT calculations is the primary objective and result of this method.

In the computation of the volume dependence of $C_{q,0}$ and $\eta$, available experimental data for atomic positions and lattice constants at various $T,V$ are used to prepare a set of ”frozen” crystal cells, each with lattice vectors and equilibrium atomic positions intended to represent the system over a range of temperatures, as shown in figure 1 . In each of the prepared cells, the atomic positions are optimized and the EFG tensor components are computed using the GIPAW method in the optimized structure, yielding the full EFG tensor at each volume, the principle axes of which comprise the ”space-fixed” coordinate system, against which the dynamical system’s rotating EFG principle axes are compared in order to ascertain $\kappa_{T}$ as in equation (12).

Inclusion of the effect of lattice dynamics requires an MD simulation in each of the prepared cells, from which the evolution of the EFG axes in a system undergoing internal motion is computed. The MD simulation generates a set of atomic configurations in each of the optimized structures, intended to represent the system at the temperature which, according to experimental data, corresponds to the cell’s volume at atmospheric pressure. The EFG tensors and the principle axes are computed in each of the perturbed structures generated by MD, and the rotation angles $\theta_{X}$, $\theta_{Y}$ which map the directions of the space-fixed axes to the updated axes are extracted.

Because intrinsic rotations do not commute we will, without loss of generality, transform the axes in the order $X^{\prime}-Y-Z$ to compute $\theta_{i}$. The average of the squares of each $\theta_{x}$ and $\theta_{y}$ over the MD trajectory determine the dynamic normalization factor $\kappa_{T}=1-\frac{3}{2}\left(\langle\theta_{X}^{2}\rangle+\langle\theta_{Y}^{2}\rangle\right)-\frac{1}{2}\eta\left(\langle\theta_{X}^{2}\rangle-\langle\theta_{Y}^{2}\rangle\right)$ which corrects the space-fixed $C_{q,0}$ for non-zero temperature. As long as the $V,T$ pairs of the previous step be chosen to correspond to the experimentally determined crystal cell and atomic positions one would find in this temperature range at atmospheric pressure, the computed frequencies would represent the $T$-dependent, laboratory-observed NQR frequencies.

IV Calculation Details

All DFT calculations were carried out using the Quantum ESPRESSO software package with PBE psuedopotentials. The geometry optimizations and the NVE molecular dynamics (MD) simulations were carried out with a cutoff energy of 60 Ry, and the EFG calculations were found to require a higher cutoff energy of 100 Ry. For both the MD and EFG calculations, an SCF-convergence threshold of 10*-7* Ry was used. The MD integration time-step was made as long as possible in both systems by constraining some bond lengths and bond angles in each system. In $\alpha$-paradichlorobenzene ($\alpha$-pCl${}_{2}\phi$), 40 constraints were applied, leaving the chlorine’s degrees of freedom as well as the translational and rotational degrees of freedom of the rigid benzene ring structure unconstrained. In HMX 16 constraints were applied, freezing only the H-C bond-lengths. The MD time-step for $\alpha$-pCl${}_{2}\phi$ was 120 a.u. (5.8 fs) and for HMX 80 a.u. (3.9 fs). Because chlorine is heavier than nitrogen, a longer simulation time was needed for $\alpha$-pCl${}_{2}\phi$ than for HMX. In both systems, the MD simulation was run for a total of 900 time steps, yielding 900 configurations used for the averaging of the EFG tensor in all cases. In both systems, a Monkhorst-Pack K-mesh was used, 2x4x4 in the case of $\alpha$-pCl${}_{2}\phi$, and 4x2x4 in the case of $\beta$-HMX.

V Structure Selection

The present method is tested for two monoclinic molecular crystals $\alpha$-pCl${}_{2}\phi$ and $\beta$-HMX. These systems were chosen because of their relative structural simplicity as well as the great body of existing research on these systems, systems, to which the calculations could be compared for the purpose of validating the methods. Each system is relatively simple, having one or few inequivalent nuclear sites and only two molecules in the unit cell. Both systems contain relatively few atomic species and a small total number of atoms in the cell: There are only 24 atoms in the unit cell of $\alpha$-pCl${}_{2}\phi$ and the 56 atom $\beta$-HMX structure was the most tractable of the popular explosives (TNT, for example, contains 168 atoms in the cell). Structural simplicity reduces computation time, a convenience for convergence testing for the many DFT calculations required by the current method. Performing the study on two quadrupolar nuclei 35Cl and 14N was also a motivation to choose one system containing each specie 35Cl and 14N.

Both systems have a long and extensive history of experimental results to which these calculations can be compared. $\alpha$-pCl${}_{2}\phi$ is the subject of some of the very earliest NQR research in history Grechishkin (1959) including on the temperature dependence of NQR Vanier (1960) and was the subject of the seminal Bayer paper Bayer (1951). The volume and stress dependence of its NQR frequencies has been characterized Kushida et al. (1956) Zamar and Brunetti (1988), and has persisted through time as a standard sample for testing new NQR techniques Sullivan (1971). The system’s intermolecular interactions and dynamics, atomic positions in the crystal and dependence of the lattice vectors on temperature have been well characterized. Frasson et al. (1959) Wheeler and Colson (1976) Estop et al. (1997)

NQR frequencies in and internal motions in $\beta$-HMX have been studied extensively Landers et al. (1981) Marino et al. (1982) Karpowicz and Brill (1983) Karpowicz and Brill (1984) Buess and Caulder (2004). As an explosive material of great interest, there is also an abundance of existing computational work associated with HMX including simulations of its internal motions Allis et al. (2006) Sewell et al. (2003). and recent DFT calculations of vibrational properties of $\beta$-HMX which have characterized the efficacy of local density approximation (LDA) vs the generalized gradient approximation (GGA) in the system. Wu et al. (2011) The $T$-dependence of the lattice vectors has been determinedDeschamps et al. (2011) and atomic positional data across a range of temperatures is available from the Cambridge Structural Database (CSD) in CIF form. Groom et al. (2016)

V.1 $\alpha$-paradichlorobenzene

Paradichlorobenzene has three polymorphs at room temperature. In this study, we examine the alpha phase. The cell is monoclinic, and contains two crystallographically equivalent molecules. Each molecule contains two equivalent 35Cl sites, making only a single inequivalent 35Cl nucleus in this structure. 35Cl has spin $\frac{3}{2}$, for which there is only one transition regardless of the magnitude of $\eta$. As a result, only a single NQR line is expected in this structure. 222 We disregard the additional line one expects due to the heavy isotope 37Cl, which we did not account for here despite its natural abundance of 24.4%. The NQR frequencies of the two isotopes should have the ratio of their respective quadrupole moments, $\nu_{{}^{35}\text{Cl}}/\nu_{{}^{37}\text{Cl}}=Q_{{}^{35}\text{Cl}}/Q_{{}^{37}\text{Cl}}=1.2704$, under normal conditions where the nuclei are not subject to significant deformations. Beginning with the structure determined by Estop et. al Estop et al. (1997), which was specified at 293K, crystal cells corresponding to the known cell parameters at five lower temperatures were generated. To obtain the atomic positions at other temperatures, the experimental cell parameters Wheeler and Colson (1976) were used to transform the atomic coordinates into the smaller cells at lower temperature. This has the undesired effect of shrinking the bond lengths, so a full geometry optimization was performed after each incremental scaling.

V.2 $\beta$-HMX

$\beta$-HMX is a monoclinic polymorph of HMX that is used as an explosives in military applications. HMX has a a pseudo-octagonal ring structure with NO2 substituents. The structure has two crystallographically equivalent molecules in the cell. Each molecule contains eight nitrogen sites, four of which are inequivalent. The four nitrogens in each molecule belonging to NO2 groups have very low NQR frequencies Landers et al. (1981). The amine nitrogens, however have NQR frequencies above 5 MHz and are therefore of greater practical interest. There are two inequivalent amine nitrogens in each molecule, differing most obviously by whether the NO2 substituent is axial or equatorial, which are descriptors of the angle of inclination the N=N bond makes with the molecular pseudoplane. Because 14N has spin 1, there are three NQR frequencies at each site in the case of an axially asymmetric field gradients. The asymmetry parameters are found to be substantial (around 0.5) for both the axial and equatorial amine nitrogens.

VI Results

The constant volume KBB model suggests a fitting function for $T$-dependence of the NQR frequencies of the form

[TABLE]

where $\nu_{0}=a_{0}$ represents the static lattice NQR frequency and the other terms result from the internal motions which give rise to equation (II.4). However, because this model is a constant volume theory, it has long been understood Wang (1955) Kushida et al. (1956) Das and Hahn (1958) the model may not suffice to describe the $T$-dependence of NQR frequencies outside of narrow temperature ranges because the static lattice coupling constant may vary significantly with cell volume. Furthermore, if the thermal expansion is anisotropic, as is the case in both systems under study in this work, variation of the cell parameters over temperature can impact $\eta$, resulting in a complicated implicit T-dependence of the parameter $\nu_{0}=a$ mediated by the thermal expansion. In the case of $\beta$-HMX, the static lattice coupling constant was not found to depend strongly on the cell volume, and consequently, the standard KBB fitting function above was used for the fits in the case of HMX. However, calculations presented in figure 5 show a strong linear dependence of $C_{q,0}$ on $V$ in $\alpha$-pCl${}_{2}\phi$ that is not present in $\beta$-HMX, whose thermal expansion is not as significant to the local EFGs. Examination of table LABEL:table:cell-parameters shows $V$ is linear in $T$ and calculations show that $\nu_{0}$ is linear in $V$; therefore, the supposition that $\nu_{0}$ is linear in T is justified and letting $\nu_{0}=a_{1}^{\prime\prime}T$ implies the appearance of a quadratic term $a_{2}T^{2}$ in equation (14)

[TABLE]

Comparison of the absolute frequency predictions, as well as the fitting parameters arising from the experimental data and the calculations may help evaluate the success of this method in both of these structures.

VI.1 $\beta$-HMX

The calculations systematically overestimate the NQR frequency in both axial and equatorial sites, likely due to error in the determination of $C_{q,0}$ since the magnitude of $\phi_{ZZ}$, defined as the largest eigenvalue of the EFG tensor, is highly sensitive to the DFT simulation parameters chosen. In particular, the suitability and quality of the pseudopotentials utilized to compute the charge densities were found to significantly impact the computed EFG tensors, even when the cutoff energies and other simulation parameters are converged independently. The degree of absolute error is reflected directly by the fitting term $a_{0}$. shown in table 3. In a static lattice, small changes in the atomic positions can also affect $C_{q,0}$ significantly; however, in this case, error introduced by atomic positional discrepancy is reduced by averaging $C_{q,0}$ over the MD trajectory. In principle one could obtain $C_{q,0}$ by performing a very well converged geometry optimization and computing the EFGs in the output. The equatorial calculations overestimate the NQR frequencies by twice as much as the axial calculations, yielding error percentages in the constant fitting parameter $a_{0}$ at 4%, and 2% for equatorial and axial sites, respectively. Consequently, the calculations underestimate the absolute difference between the axial and equatorial NQR frequencies by about a factor of $\frac{1}{2}$. The electronic environments of the two sites types differ primarily due to the orientation of the substituent NO2 molecule with respect to the molecular pseudo-plane, to which the DFT is sensitive. Using the fitting function $\nu\left(T\right)=a_{0}+a_{1}T+a_{-1}/T$ the experimental data from Landers et al. (1981) for showed the axial amines to have generally NQR frequencies about 200kHz higher than the equatorial sites, and a slightly sharper fall-off with temperature. No additional quadratic term was used to fit either experimental or calculated frequencies because the computed static lattice values of $C_{q,0}$ did not strictly increase with increasing volume, as shown in figure 5.

VI.2 $\alpha$-paradichlorobenzene

In this structure, the calculations systematically underestimate the NQR frequency by about 500KHz, as shown in figure 7. The cause of the underestimation is again likely error in the estimation of $C_{q,0}$ which has a large impact on the absolute errors in the calculation. The T-dependence of NQR frequencies $\alpha$-paradichlorobenzene was fit using a modified KBB fitting function given by equation (15), $\nu\left(T\right)=a_{0}+a_{1}T+a_{2}T^{2}$. The experimental NQR data resulting in fitting parameters $a_{0}=34.8437$ MHz, $a_{1}=-5.2424\times 10^{-4}\text{ MHz}\cdot T^{-1}$, and $a_{2}=-4.6970\times 10^{-6}\text{ MHz}\cdot T^{-2}$. The corresponding fitting parameters $a_{i}$ from the fits of the calculated NQR frequencies are $a_{0}=34.2932$ MHz, $a_{1}=-5.5559\times 10^{-4}\text{ MHz}\cdot T^{-1}$, and $a_{2}=-4.2894\times 10^{-6}\text{ MHz}\cdot T^{-2}$.

The fitting parameters are summarized in table 3. The constant terms agree to within 1.6% and the dominant quadratic terms agree within 12%, while the linear terms differ by 26%. The absolute error of the NQR frequency prediction across the 100K-293K temperature range was at most 580 KHz and at its lowest 450KHz, while the amount by which the measured NQR frequency drifts over this temperature range is similarly about 450 KHz.

VII Discussion

VII.1 Advantages and Disadvantages

The method proposed has the advantage that it can be applied in a straightforward fashion to any system in which data for atomic positions and the thermal expansion of the lattice vectors are available, without the need for additional postulates or supporting experimental data such as IR or Raman studies. A disadvantage of the current application of the method is that the computational cost can be high; in the present study, an MD simulation of 900 steps followed by 900 SCF/EFG calculations were used to generate a single data point on the plots in figures 6 and 7. Six and sixteen data points were computed in $\alpha$-pCl${}_{2}\phi$ and $\beta$-HMX, respectively.

VII.2 Considerations for Proper Implementation

The technique presented here is structure agnostic because the heavy-lifting is done by the DFT calculations. In principle, the results will be accurate in other systems provided the MD simulations faithfully represent the structure’s dynamics at the desired temperatures, and the subsequent sets of EFG calculations are carried out with sufficient simulation parameters for accuracy, requiring some care in the application of the simulations. Attention to detail regarding both the underlying physics and best practices in computing is essential for a successful implementation, especially for the management of the output data, having both considerable size and complexity.

For the EFG calculations, careful choice of pseudopotentials and the simulation parameters of the DFT calculations, such as cutoff energies and SCF convergence thresholds is needed in order to ensure convergence of the EFGs. It is also worth considering the precision of the known quadrupole moments in the nuclei under study, and recognition that the values of $Q$ used in the EFG calculations are constant here, though deformation of the nucleus by the surroundings can occur and effectively change the instantaneous value of $Q$ for the nucleus. 333NQR measurements involve direct observation of the coupling constants $eQ\phi_{ZZ}/h$, whereas in the present study, we leave $Q$ a constant and account for dynamical effects on $\phi_{ZZ}$. Values of $Q$ in standard tables are of the intrinsic quadrupole moments. But the nuclear environment may cause deformations of the nucleus that induce contributions to the quadrupole moment. Finally, One must perform the summary operation of averaging the EFG tensors in the correct manner to compute the correction to the NQR frequency computed in the static lattice.

The data generated may have considerable size on disk making computation of the $\langle\theta_{i}\rangle$ cumbersome without some consideration taken in the efficient storage and processing of the default simulation outputs. On occasion, the convention $|\phi_{ZZ}|\geq|\phi_{YY}|\geq|\phi_{XX}|$ will cause a sudden swap in the names EFG’s eigenvectors if the magnitude of $\phi_{XX}$ momentarily overtakes $\phi_{YY}$, leading to a discontinuity in the time series of the $\theta_{i}$. This can be remedied by checking that the moving principal axes coordinate system remains right handed.

VIII Conclusion

The calculations underestimate the NQR frequencies in $\alpha$-pCl${}_{2}\phi$ and overestimate them in $\beta$-HMX. The overestimation is about twice are large in the equatorial sites than for the axial ones. The calculated and experimental fits agreed reasonably well when the experimental data was restricted to the same temperature range covered by the calculations. In $\alpha$-pCl${}_{2}\phi$ the fitting of the calculated NQR frequencies match the experimental fits to 1.56%, and 12.9% for the constant $a_{0}$ and quadratic $a_{2}$ terms respectively, but the linear term $a_{1}$ displayed error in excess of 25%.

In $\beta$-HMX the experimental and calculated constants $a_{0}$ agreed to 2.09% and 4.00% for axial and equatorial sites, respectively. The errors in the linear terms were higher at 14.3% and 19.9%. The $1/T$ fitting terms $a_{-1}$ agreed to less than 0.5% in the equatorial case but the error was 22.3% in the axial case. The $a_{-1}$ term likely requires data at lower temperatures to be fit reliably, as the lowest temperature calculated in $\beta$-HMX was 123K.

The method presented may be applied to arbitrary solids but more work is needed to extend the method to larger systems in a scalable manner. Repeated application of the general method is needed to confirm its robustness, for which automation of both the input data preparation and output data processing would be of great use since the datasets generated can reach an unruly size on disk in systems of even modest size, without clever implementation. On-the-fly data pruning, parsing and processing will likely be needed to extend this method to very large systems.

The DFT computations carried out in this procedure produce a rich dataset that may be exploited in various ways that may lend improved results or computational efficiency, such as in the manner of tensor averaging used to compute $\kappa_{T}$. There is also the potential to train statistical models to map atomic positions to EFGs which, if effective, would substantially increase efficiency as the molecular dynamics simulations are significantly cheaper than the SCF and EFG computations in terms of CPU time.

IX Acknowledgements

This work is supported in part by the National Science Foundation grant NSF-DMR-1157490 and Department of Energy grants DE-FG02-02ER45995. The authors extend gratitude to NERSC for computational resources and its support staff for assistance with compilations and automation procedures.