Theoretical predictions of the two-dimensional solid-state NMR spectra: A case study of the 13C–1H correlations in metergoline

Chemical Physics Letters 586 (2013) 56–60

Contents lists available at ScienceDirect

Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett

Theoretical predictions of the two-dimensional solid-state NMR spectra: A case study of the 13C–1H correlations in metergoline Jirˇí Czernek ⇑, J. Brus Institute of Macromolecular Chemistry, Academy of Sciences of the Czech Republic, Heyrovsky Square 2, 162 06 Praha 6, Czech Republic

a r t i c l e

i n f o

Article history: Received 27 June 2013 In final form 4 September 2013 Available online 20 September 2013

a b s t r a c t A new method for the treatment of data from multidimensional solid-state NMR investigations is described. It approximates the theoretical NMR chemical shifts from the chemical shielding values obtained by first-principles calculations and subsequently treats these results to quantify the similarity between predicted and experimental chemical shift correlations. The test case of this approach is performed for the measured and several sets of computed 13C–1H heteronuclear correlations in the polymorphic form I of metergoline, which is relatively large, pharmaceutically active system. The proposed protocol is general, however, and it can be immediately applied to study other compounds and nuclei. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction Modern methods of the solid-state nuclear magnetic resonance (SSNMR) can be usefully employed to study a wide range of problems in chemistry, structural biology and material science [1]. Many of the SSNMR contributions to these fields have involved the measurements of the two-dimensional (2D) spectra (see Ref. [2] for a comprehensive review, which describes the underlying pulse sequences together with their applications). Probably the most common 2D SSNMR experiments, also discussed in detail in Ref. [2], correlate the 13C and 1H chemical shifts. For a system containing these two atoms it is usually necessary to carry out the 13 C–1H correlation experiments to obtain the assignment of the SSNMR signals to specific nuclei in the investigated solid-phase structure, which is a prerequisite for virtually any subsequent investigations. As pioneered by the groups of Harris [3] and Emsley [4] (see also the most recent reviews [5,6]), the process of assigning the resonances in the measured spectra can benefit from an application of the first-principles calculations of the NMR chemical shielding property [1] of the known or tentative periodic arrangements of solids. These calculations were also used to predict the information obtainable from measurements of the 13C–1H correlations. The most notable example involved the gauge-including projector augmented-wave (GIPAW) [7,8] computations for campho[2,3-c]pyrazole, which were indispensable for the assignment of the six independent molecules in its unit cell [9]. However, these and other available [10,11] theoretical predictions of SSNMR heteronuclear correlations were not quantitatively evaluated, but only visually compared to the positions of cross-peaks in measured spectra. In this Letter, a simple but general measure of the quality ⇑ Corresponding author. Fax: +420 296809410. E-mail address: [email protected] (J. Czernek). 0009-2614/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cplett.2013.09.015

of predicted heteronuclear correlations is proposed. Its form is derived from a thorough analysis of the combined SSNMR experimental and density-functional theory (DFT) study of the 13C–1H chemical shift correlations in metergoline, a medium-sized organic molecule and potent drug [12]. This measure, which is defined in Section 2.3, is used to investigate three aspects of the reproducibility of the 2D SSNMR spectra by DFT-based computations. First, the influence of crystal-lattice effects upon the 13C and 1H chemical shifts is described by a careful comparison of the quality of the results predicted for various structural models of metergoline. Second, two plane-wave DFT computer codes, namely, CASTEP [13] and Quantum ESPRESSO [14] (QE) are employed and their performance in the reconstruction of the 13C–1H correlation spectra is critically compared. Third, by applying the QE program and the van der Waals density-fitting (vdW-DF) approach of de Gironcoli and coworkers [15], DFT calculations are performed with the full inclusion of London dispersion force, whose role in simulations of 2D SSNMR spectra is assessed. Additionally, the quality of predictions obtained using the CASTEP program and three exchangecorrelation functionals (specified in Section 2.2) is evaluated for the 13C and 1H sets of data. The confrontation of the measurements with theoretical predictions [16], described here for metergoline, is of importance for high-precision structural studies of active pharmaceutical ingredients (for recent examples, see [17] and references therein) and their cocrystals [18]. 2. Methods 2.1. Experimental Metergoline (C25H29N3O2, see Figure 1; IUPAC name: Benzyl N-{[(6aR,9S,10aR)-4,7-dimethyl-6,6a,8,9,10,10a-hexahydroindolo[4, 3-fg]quinolone-9-yl)methyl}carbamate; CAS number: 17692-51-2)

J. Czernek, J. Brus / Chemical Physics Letters 586 (2013) 56–60

57

Figure 1. The atom numbering scheme adopted for an evaluation of the NMR data of metergoline. Exclamation marks denote experimentally unresolved sites (see Section 2.3), while ‘a’ and ‘e’ are used to discriminate between protons in the axial and equatorial position, respectively.

was used as obtained from TEVA Industries s.r.o., Czech Republic. It was ascertained by the X-ray powder diffraction (XRPD) analysis that the sample exclusively contained the polymorphic form I [19,20]. The SSNMR experiments were run using 4 mm double-resonance probe head of a Bruker Avance 500 WB/US spectrometer (the carrier frequency of 500.18 and 125.78 MHz for 1H and 13C nuclei, respectively). In brief, following the standard variable-amplitude cross-polarization measurement of the 13C signals, performed at the magic-angle spinning frequency of 12 kHz, a series of 2D experiments was carried out. The cross-polarization polarization-inversion experiment [21] applied with the inversion period of 55 ls enabled the discrimination of methylic, methylenic, primary and quaternary carbons of metergoline. This was followed by the frequency-switched Lee–Goldburg measurements [22], which were performed with three mixing periods (0.1, 0.3 and 3.0 ms) thus allowing to distinguish the spatial proximities of the respective 13C–1H spin-pairs [23]. Additional experimental details can be found in Ref. [20]. 2.2. Computational Both the CASTEP and QE programs adopt the Kohn–Sham DFT and treat the crystalline structure as an infinite system using periodic boundary conditions, cf. Refs. [13,14,24,25] for details. These two programs were employed to carry out the full optimizations of internal coordinates of the polymorphic form I of metergoline [20]. The unit cell parameters were fixed at values found in the single-crystal X-ray diffraction study [19], and the initial coordinate set was taken also from Ref. [19], in all five geometrical optimizations specified below. Thus, three generalized gradient approximation (GGA) DFT exchange-correlation functionals, i.e., PBE [26], PW91 [27] and revPBE (‘revised PBE’) [28], were applied together with the ‘Fine’ level of settings of the CASTEP 5.0 program. In particular, the 3  2  5 Monkhorst–Pack [29] (MP) grid of 15 kpoints, and the plane-wave cut-off energy of 550 eV, were adopted. During the two geometrical optimizations performed using the Quantum ESPRESSO 4.3.1 package, the convergence criteria were set to tight values of 1.0D-5 and 1.0D-6 Hartree atomic units for the total crystal-lattice energy and for its gradient, respectively. The same MP grid as in the CASTEP calculations was used, while the cut-off value was set to 1361 eV. As for the functionals employed, one was the standard PBE, while the other (the ‘revPBEvdW-DF’ recently implemented in the QE distribution by de Gironcoli and coworkers [15]) combined the revPBE exchange with the local density correlation term, which perturbatively included dispersion effects [30,31]. The average force per atom was 3.7046e04, 7.1635e-04, 5.3582e-04, 4.2271e-05 and 1.9174e-05 for the final structures obtained by the PBE (CASTEP), PW91 (CASTEP), revPBE (CASTEP), PBE (QE) and revPBE-vdW-DF (QE) optimizations, respectively. The predictions of the NMR chemical shielding tensors were performed using the same programs as those employed for obtain-

Figure 2. The dimer model with the deshieldied axial-H4 atom shown as the white ball. See the text for details.

ing the optimized coordinates, i.e., the NMR-CASTEP module (the ‘Fine’ settings) and the QE-GIPAW version 4.3.1 (the default thresholds). In the case of GGA functionals, the GIPAW approach [7,8] was used. In the QE suite of codes, the revPBE-vdW-DF calculations are only implemented within the projector augmented-wave (PAW) formalism [32], and they require an application of the pseudopotential generation routine, which is available from the PSLibrary project [33]. The pseudopotentials thus obtained for C, H, N and O can be received from the corresponding author upon request. The supramolecular cluster [34,35] was built to explicitly model the influence of intermolecular interactions on the 1H and 13C chemical shielding of the polymorphic form I of metergoline. The final coordinate set provided by the PBE geometrical optimization in CASTEP was used and two neighboring molecules were clipped out, which were positioned in the direction of the lattice vector b of the metergoline crystal [19]. This structure is referred to as ‘dimer’ and shown in Figure 2. In addition, an isolated molecule (‘monomer’) was considered. The difference between the dimer and monomer values of the chemical shifts can be usefully employed for an assessment of the role of intermolecular interactions upon NMR parameters [36,37]. For these two structures, the chemical shielding tensors were calculated using the GAUSSIAN 09 package [38] by adopting the GIAO [39] method to overcome the gauge problem and by applying the B3LYP [40,41] combination of DFT exchange-correlation functionals together with the TZ2P [42] basis set. 2.3. Data evaluation From values of the principal elements, rxx, ryy, rzz, of the chemical shielding tensors, computed by the methods specified above, the absolute isotropic NMR chemical shielding, r, r = (rxx + ryy + rzz)/3, was obtained. The r results found for the 1H and 13C sites of metergoline were compared to the measured values of the isotropic NMR chemical shifts, d, without any attempts on the referencing of calculated data [3]. As the experiments could not distinguish the signals of the protons in ortho positions of the

58

J. Czernek, J. Brus / Chemical Physics Letters 586 (2013) 56–60

phenyl ring, and of the N1- and N6-methyl protons, the corresponding r values were averaged, thus reducing the number of investigated 1H data-points from 29 to 24 (see Figure 1). Analogously, because the two carbons in meta positions were not experimentally resolved, the reduction was also performed for the 13C results (from 25 to 24). The reduced sets of r were linked to their measured counterparts d by finding the coefficients p and q from the least-squares fit to the functional form:

r ¼ p  d þ q:

ð1Þ

The bivariate data, representing the 13 carbons with 18 hydrogens attached to them, were considered in modeling the 13C–1H correlation spectra. Thus, first the coefficients a, b, c and d are obtained, which minimize:

mindðCÞ;a;b

X13 i¼1

2

ða  dðCÞi þ b  rðCÞi Þ

ð2Þ

and

mindðHÞ;c;d

X18 j¼1

2

ðc  dðHÞj þ d  rðHÞj Þ :

ð3Þ

where r(C)i and d(C)i denote, respectively, the 13C chemical shielding and 13C chemical shift of the carbon atom i, and r(H)j and d(H)j accordingly stand for the 1H chemical shielding and 1H chemical shift of the hydrogen atom j. Then the theoretical 13C chemical shift of the carbon atom i, e(C)i, and the theoretical 1H chemical shift of the hydrogen atom j, e(H)j, are extracted from:

eðCÞi ¼ a  rðCÞi þ b

ð4Þ

and

eðHÞi ¼ c  rðHÞi þ d:

ð5Þ

It is obvious that the experimental 13C–1H correlation spectrum of metergoline can be described by specifying the 18 pairs of chemical shift values [d(C)i; d(H)j]. Their predicted counterparts [e(C)i; e(H)j] can be employed to estimate an agreement between theory and measurement as follows. For n such pairs and for each j = 1, 2, . . ., n, the differences are calculated:

qðHÞj ¼ dðHÞj  eðHÞj

ð6Þ

and

pðCÞj ¼ dðCÞj  eðCÞj

ð7Þ

where indexing the carbons with j designates the direct correspondence between the hydrogen(s) bound to a given carbon. Thus, two Table 1 The statistical data describing an agreement between the predicted

13

identical values of p(C)j are to be used for the carbon with two protons attached to it (in the case of metergoline, this concerns C4, C7, C9, C17 and C21 atoms, cf. Figure 1), while all q(H)j are generally distinct. Then the covariance between these differences, sqp, is obtained as:

P Sqp ¼

P j

qðHÞj pðCÞj 

j

qðHÞj

P j

pðCÞj



n

ð8Þ

n1

and quantifies the similarity of the original sets [d(C)i; d(H)j] and [e(C)i; e(H)j] (and, by extension, [r(C)i; r(H)j]). 3. Results and discussion 3.1. Theory and experiment for the

13

C and 1H NMR chemical shifts

As described above, for 24 carbons and 24 protons of the polymorphic form I of metergoline the chemical shifts were unambiguously assigned, and the chemical shieldings computed by various theoretical approaches (all these values are available from Supplementary information files C.txt and H.txt). Using the linear fits expressed by Eq. (1), the calculated data were confronted with their experimental counterparts. The results of this evaluation are summarized in Table 1. They clearly show an excellent agreement between the predictions obtained from the five methods, which were employed to properly treat the solid-phase environment, and the measured chemical shifts. Specifically, the average values of one standard deviation of the fits of these theoretical results to experiment amount to 1.02 and 0.22 ppm for 13C and 1H data, respectively, and in the case of the maximum absolute deviation they accordingly equal to 2.52 and 0.44 ppm. Consistent with findings reported recently for the 13C chemical shielding in peptides [43] and in thiamine-based systems [44], the results calculated in CASTEP with the PBE, PW91 and revPBE functionals are all in agreement, including the respective values of the slope (p) and intercept (q) in Eq. (1). Expectedly, the predictions obtained for an isolated molecule clipped out from the fully optimized crystal of metergoline are dubious, especially for the shielding of protons (see Table 1). The most drastic example is the overestimation by more than two ppm of the chemical shielding of the axial-H4 nucleus, which exhibits a small value of the chemical shift (d = 0.27 ppm) due to the secondary structural effects upon the chemical shielding parameters [45,46]. An application of the simple cluster model, graphically presented in Figure 2, only qualita-

C (in parentheses, 1H) NMR isotropic chemical shieldings and measured chemical shifts.

Chemical shielding calculation

Geometry

Slope

Standard error of slope

Intercept (ppm)

Standard error of intercept (ppm)

Standard deviation (ppm)

Average absolute deviation (ppm)

Maximum absolute deviation (ppm)

Adjusted R2

B3LYP-GIAO/TZ2P (GAUSSIAN 09) for the monomer B3LYP-GIAO/TZ2P (GAUSSIAN 09) for the dimer PBE (CASTEP) for the crystal PW91 (CASTEP) for the crystal revPBE (CASTEP) for the crystal PBE (Quantum Espresso) for the crystal revPBE-vdW-DF (Quantum Espresso) for the crystal

Taken from the crystal (PBE in CASTEP)

1.0412 (0.8518)

0.0111 (0.0643)

175.20 (31.17)

1.13 (0.35)

2.26 (0.76)

1.74 (0.57)

4.56 (2.15)

0.99737 (0.88366)

Taken from the crystal (PBE in CASTEP)

1.0440 (0.9813)

0.0103 (0.0479)

175.58 (32.07)

1.04 (0.26)

2.08 (0.57)

1.61 (0.38)

4.15 (1.39)

0.99778 (0.94787)

Optimized: PBE in CASTEP Optimized: PW91 in CASTEP Optimized: revPBE in CASTEP Optimized: PBE in Quantum Espresso Optimized: revPBEvdW-DF in Quantum Espresso

1.0292 (1.0456) 1.0273 (1.0446) 1.0187 (1.0437) 1.0419 (1.1016) 1.0102 (1.0533)

0.0047 (0.0184) 0.0048 (0.0190) 0.0052 (0.0191) 0.0048 (0.0168) 0.0057 (0.0187)

172.52 (31.99) 172.20 (32.04) 173.35 (32.18) 166.59 (25.09) 159.63 (25.39)

0.48 (0.10)

0.96 (0.22) 0.95 (0.22) 1.06 (0.22) 0.97 (0.20) 1.16 (0.22)

0.66 (0.18)

2.35 (0.42)

0.64 (0.19)

2.46 (0.45)

0.77 (0.19)

2.46 (0.44)

0.66 (0.15)

2.70 (0.44)

0.92 (0.18)

2.66 (0.44)

0.99951 (0.99295) 0.99952 (0.99243) 0.99940 (0.99238) 0.99952 (0.99466) 0.99926 (0.99279)

0.48 (0.10) 0.53 (0.10) 0.48 (0.09) 0.58 (0.10)

59

J. Czernek, J. Brus / Chemical Physics Letters 586 (2013) 56–60

Table 2 The level of agreement between the measured 13C NMR isotropic chemical shifts of 13 carbon nuclei and the corresponding isotropic chemical shieldings predicted by various methods (in parentheses, analogous results for 18hydrogen atoms), and the values of the covariance defined in the text. Chemical shielding calculation

Geometry

Slope

Intercept (ppm)

Standard deviation (ppm)

Adjusted R2

Covariance

B3LYP-GIAO/TZ2P (GAUSSIAN 09) for the monomer PBE (CASTEP) for the crystal

Taken from the crystal (PBE in CASTEP)

revPBE-vdW-DF (Quantum Espresso) for the crystal

Optimized: revPBE-vdW-DF in Quantum Espresso

0.99844 (0.84972) 0.99940 (0.99229) 0.99935 (0.99367) 0.99893 (0.99062)

0.6476

Optimized: PBE in Quantum Espresso

169.67 (32.54) 168.00 (30.18) 160.56 (22.52) 159.15 (23.89)

1.47 (0.93)

PBE (Quantum Espresso) for the crystal

0.9815 (1.0272) 0.9743 (0.9409) 0.9657 (0.8952) 1.0011 (–0.9390)

Optimized: PBE in CASTEP

tively improves the results (cf. Table 1). For example, in spite of the inclusion of a neighboring molecule lying about 3.0 Å atop of the axial-H4 nucleus, an error in its predicted chemical shielding still exceeds 1 ppm. Nevertheless, in order to test thoroughly the method proposed for the simulation of the 13C–1H correlation spectra, also the results from the calculation describing an isolated metergoline molecule were utilized. Thus, they were included in, formally speaking, the experiment-to-theory fits, described by Eqs. (2) and (3), together with the values obtained from several highquality methods (this procedure serves the purpose of a straightforward extraction of the theoretical chemical shifts from Eqs. (4) and (5)). Table 2 summarizes the key parameters of these linear correlations, which were obtained for the data describing carbons C2, C4, C5, C7, C8, C9, C17, C21 and C25, and for the 18 protons attached to them (see Figure 1). With the anticipated exception of the results for the ‘monomer’ protons, an agreement between the measured and calculated sets is satisfactory. Notably, the goodness-of-fit statistics for this analysis provides highly similar values to those obtained from the model described by Eq. (1) (cf. Tables 1 and 2). As a consequence, the quality of results for the 13 carbons specified above, and for their adjacent protons, is not markedly different from that found for the full set of 13C and 1H nuclei. Hence, in the case of metergoline molecule, the differences given by Eqs. (6) and (7) should provide reasonable estimates of inaccuracies in the predicted 2D spectra. 3.2. Modeling the

13

0.90 (0.21) 0.95 (0.19) 1.22 (0.23)

0.0692 0.0413 0.0693

Figure 3. The 13C–1H chemical shift correlations in metergoline as obtained from the experiment (open circles), the monomer model (full circles), and the GIPAW PBE calculations in CASTEP (asterisks), together with the corresponding values of the covariance. See the text for details.

C–1H correlations

Table 2 also shows the values of the covariance (cf. Eq. (8)) between the differences q(H) and p(C) of experimental chemical shifts and their theoretical counterparts, which were obtained by the four theoretical approaches with vastly differing quality of predicted heteronuclear correlations. Figures 3 and 4 present the results from selected combinations of computational methods together with the experimental data defining the 13C–1H correlations, i.e., the [d(C)i; d(H)j] pairs. As it could be inferred from a visual inspection of Figure 3, the peak positions are poorly predicted if crystal-lattice effects are included only indirectly, i.e., by the B3LYP-GIAO/TZ2P calculations on an isolated molecule taken from the optimized crystal structure. Importantly, the resulting value of the covariance, 0.65, is one order of magnitude higher than for the data obtained by the GIPAW and PAW approaches (cf. Table 2). In order to see what effects the respective methods have on the predicted 2D correlation spectra, attempts can also be made to compare much smaller differences in the covariance. Thus, it appears that the calculations carried out in the QE package using the PBE functional were the most successful from those tested here, with the covariance as low as 0.04. However, it should be kept in mind that the convergence criteria were tightened with respect to the default values for geometry optimization (see Section 2.2), and the covariance resulting from analo-

Figure 4. The 13C–1H chemical shift correlations in metergoline as obtained from the experiment (full circles) and from the Quantum ESPRESSO calculations using the GIPAW PBE (open circles) and the PAW revPBE-vdW-DF (asterisks) approaches, together with the corresponding values of the covariance. See the text for details.

60

J. Czernek, J. Brus / Chemical Physics Letters 586 (2013) 56–60

gous computations performed in CASTEP is only slightly increased, so the quality of these two predictions, while being excellent, can be considered the same. As for there revPBE-vdW-DF strategy, it was successfully applied to systems with strong intermolecular interactions, in particular, cholesterol crystals [47] and polymorphs of glycine under pressure [48]. Taking into account an absence of hydrogen bonds in the polymorphic form I of metergoline [20], it is not surprising that the 13C–1H correlations predicted by this novel approach do not significantly differ from their PBE counterparts (see Table 2 and Figure 4). Of course, it might be advantageous to apply the revPBE-vdW-DF method (or its extensions [49]) in simulations of some other multidimensional SSNMR spectra. 4. Conclusion First-principles calculations of the SSNMR parameters are becoming an important component of experimental investigations, and the possibilities of their applications were extended by developing a simple strategy for the quantification of similarity between heteronuclear correlations of the chemical shifts. The test case of this approach was performed for the set of measured and predicted 13 C–1H correlations in a relatively large system, i.e., the polymorphic form I of metergoline. Employing these data, the quality of spectral predictions from several methods of an inclusion of the intermolecular effects upon the chemical shifts was evaluated. Importantly, the PBE geometrical optimization of the atomic positions within the crystal followed by the GIPAW chemical shielding prediction was found to lead to the results, which would have enabled the correct (and independent) assignment of the investigated signals. Acknowledgements Computational resources were partially provided under the programme ‘Projects of Large Infrastructure for Research, Development, and Innovations’ LM2010005, and in the Center CERIT Scientific Cloud, part of the Operational Program Research and Development for Innovations, reg. no. CZ. 1.05/3.2.00/08.0144. We are grateful to Dr. Michal Hušák, Institute of Chemical Technology in Prague, for the XRPD analysis of metergoline and for the technical assistance with CASTEP calculations, and to Dr. E. Küçükbenli, École Polytechnique Fédérale de Lausanne, for detailed instructionson the usage of the PSlibrary.0.3.0. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.cplett.2013.09. 015. References [1] M.J. Duer, Introduction to Solid-State NMR Spectroscopy, Blackwell, Oxford, 2004.

[2] S.P. Brown, Solid State Nucl. Magn. 41 (2012) 1. [3] R.K. Harris, P.Y. Ghi, Magn. Reson. Chem. 44 (2006) 325. [4] N. Mifsud, B. Elena, C.J. Pickard, A. Lesage, L. Emsley, Phys. Chem. Chem. Phys. 8 (2006) 3418. [5] C. Bonhomme et al., Chem. Rev. 112 (2012) 5733. [6] T. Charpentier, Solid State Nucl. Magn. 40 (2011) 1. [7] C.J. Pickard, F. Mauri, Phys. Rev. B 63 (2001) 245101. [8] J.R. Yates, C.J. Pickard, F. Mauri, Phys. Rev. B 76 (2007) 024401. [9] A.L. Webber, L. Emsley, R.M. Claramunt, S.P. Brown, J. Phys. Chem. A 114 (2010) 10435. [10] R.K. Harris, P. Hodgkinson, T. Larsson, A. Muruganantham, I. Ymén, D.S. Yufit, V. Zorin, Cryst. Growth Des. 8 (2008) 80. [11] S. Cadars, A. Lesage, C.J. Pickard, P. Sautet, L. Emsley, J. Phys. Chem. A 113 (2009) 902. [12] J.M. Hooker, S.W. Kim, A.T. Reibel, D. Alexoff, Y. Xu, C. Shea, Med. Chem. 18 (2010) 7739. [13] S.J. Clark, M.D. Segall, C.J. Pickard, P.J. Hasnip, M.J. Probert, K. Refson, M.C. Payne, Z. Kristallogr. 220 (2005) 567. [14] P. Giannozzi et al., J. Phys. Condens. Matter 21 (2009) 395902. [15] A.R. Ferreira, E. Küçükbenli, A.A. Leitao, S. de Gironcoli, Phys. Rev. B (2011) 235119. [16] M. Profeta, F. Mauri, C.J. Pickard, J. Am. Chem. Soc. 125 (2003) 541. [17] A.S. Tatton, T.N. Pham, F.G. Vogt, D. Iuga, A.J. Edwards, S.P. Brown, CrystEngCom 14 (2012) 2654. [18] M. Khan, V. Enkelmann, G. Brunklaus, J. Am. Chem. Soc. 132 (2010) 5254. [19] E.F. Serantoni, P. Sabatino, L.R. Sanseverino, G.M. Sheldrick, Acta Crystallogr. B33 (1977) 2899. [20] M. Hušák et al., Struct. Chem. 19 (2008) 517. [21] X. Wu, T. Burns, K.W. Zilm, J. Magn. Reson. A 111 (1994) 29. [22] B.J. van Rossum, C.P. de Groot, V. Ladizhansky, S. Vega, H.J.M. de Groot, J. Am. Chem. Soc. 122 (2000) 3465. [23] J. Brus, A. Jegorov, J. Phys. Chem. A 108 (2004) 3955. [24] G. Kresse, D. Joubert, Phys. Rev. B 59 (1999) 1758. [25] M.D. Segall, P.J.D. Lindan, M.J. Probert, C.J. Pickard, P.J. Hasnip, S.J. Clark, M.C. Payne, J. Phys. Condens. Matter 14 (2002) 2717. [26] J.P. Perdew, K. Burke, Phys. Rev. Lett. 77 (1996) 3865. [27] J.P. Perdew, Y. Wang, Phys. Rev. B33 (1986) 8800. [28] B. Hammer, L.B. Hansen, J.K. Norskov, Phys. Rev. B59 (1999) 7413. [29] H.J. Monkhorst, J.D. Pack, Phys. Rev. B 13 (1976) 5188. [30] M. Dion, H. Rydberg, E. Schröder, D.C. Langreth, B.I. Lundqvist, Phys. Rev. Lett. 92 (2004) 246401. [31] G. Román-Pérez, J.M. Soler, Phys. Rev. Lett. 103 (2009) 096102. [32] P.E. Blochl, Phys. Rev. B 50 (1994) 17953. [33] A.A. Adllan, A.D. Corso, J. Phys. Condens. Matter 23 (2011) 425501. [34] C. Ochsenfeld, F. Koziol, S.P. Brown, T. Schaller, U.P. Seelbach, F.-K. Klärner, Solid State Nucl. Magn. 22 (2002) 128. [35] K. Malinˇáková, L. Novosadová, M. Lahtinen, E. Kolehmainen, J. Brus, R. Marek, J. Phys. Chem. A 114 (2010) 1985. [36] J.R. Yates, T.N. Pham, C.J. Pickard, F. Mauri, A.M. Amado, A.M. Gil, S.P. Brown, J. Am. Chem. Soc. 127 (2005) 10216. [37] J. Schmidt, A. Hoffmann, H.W. Spiess, D. Sebastiani, J. Phys. Chem. B 110 (2006) 23204. [38] M.J. Frisch et al., GAUSSIAN 09, Revision C.01, Gaussian, Inc., Wallingford, CT, 2010. [39] K. Wolinski, J.F. Hinton, P. Pulay, J. Am. Chem. Soc. 112 (1990) 8251. [40] A.D. Becke, J. Chem. Phys. 98 (1993) 5648. [41] C. Lee, W. Yang, R. Parr, Phys. Rev. B 37 (1998) 785. [42] A. Schäfer, H. Horn, R. Ahlrichs, J. Chem. Phys. 97 (1992) 2571. [43] J. Czernek, T. Pawlak, M.J. Potrzebowski, Chem. Phys. Lett. 627 (2012) 31. [44] J. Czernek, T. Pawlak, M.J. Potrzebowski, J. Brus, Chem. Phys. Lett. 555 (2013) 135. [45] D. Sitkoff, D.A. Case, Prog. Nucl. Magn. Reson. Spectrosc. 32 (1997) 165. [46] J. Czernek, J. Phys. Chem. A 107 (2003) 3952. [47] E. Küçükbenli, K. Sonkar, N. Sinha, S. de Gironcoli, J. Phys. Chem. A116 (2012) 3765. [48] R. Sabatini, E. Küçükbenli, B. Kolb, T. Thonhauser, S. de Gironcoli, J. Phys. Condens. Matter 24 (2012) 424209. [49] V.R. Cooper, Phys. Rev. B 81 (2010) 161104.