Towards quantitative prediction of proton chemical shifts in imidazolium chloride ionic liquids by computational NMR

Abstract

We present the results of a computational investigation based on molecular dynamics (MD) and density functional theory (DFT) aimed at testing the performance of different classical Force Fields (FFs) to accurately model the microscopic structure of the ionic liquids butylmethylimidazolium tetrafluoroborate and chloride and predict their NMR properties. The FFs used are a full-charge FF, a scaled-charge FF, a polarizable FF, and a recently presented virtual-site FF (Doherty, B.; Zhong, X.; Acevedo, O. J. Phys. Chem. B 2018, 122, 2962–2974). We use the NMR parameters (namely ¹H chemical shifts of the [C₄C₁im] cation) computed by DFT methods over a MD trajectory as a probe to judge the performance of the various FFs. While the agreement between experimental and calculated NMR resonances is almost quantitative for the tetrafluoroborate salt, using all four types of FFs, for the chloride salt we observe a very different performance: The full-charge and scaled-charge FFs exhibit a rather significant disagreement, while an improvement is obtained with the polarizable FF and a rather good agreement is achieved by the virtual-site FF. The high sensitivity of NMR parameters to geometrical factors allows us to pinpoint specific deficiencies of different FFs in correctly representing the hydrogen-bond between the positively charged imidazolium ring protons and the chloride anion. In turn, such analysis enables us to propose future directions for improving the performance of classical FFs for ionic liquids.

---

Introduction

Computational NMR refers to the vast array of computational protocols, mostly based on electronic density functional theory (DFT) methods, used to predict the NMR properties of putative structures of unknown compounds. Such methods were pioneered about 25 years ago by Bagno,¹ Bifulco,² Köck,³ and Sebag⁴ for organic molecules and natural substances, and were further developed later by the groups of Goodman,^5,6 and Sarotti.^7–9 Recent review papers have highlighted the state of art in this field.^10,11 In short, the NMR data (¹H and ¹³C chemical shifts and spin–spin coupling constants) of a given compound with an unknown structure are compared with the results of accurate DFT calculations of hypothetical structures, using several statistical parameters, such as the correlation coefficient of the linear fit, the mean absolute error (MAE), etc. Such a comparison makes it possible to discard or retain a given structural proposal based on the agreement between calculated and experimental values. For organic molecules and natural substances, the “structural space” to be explored is discrete, meaning that there are several covalent structures that can be envisaged from the beginning, which differ qualitatively from each other.

Recently, the same idea was proposed and applied to more complex systems such as bulk ionic liquids (ILs).^12,13 It is well known that the chemical shift of the ions constituting ILs is strongly dependent on the counterion. For instance, the ¹H resonance of the H2 ring proton of 1-butyl-3-methylimidazolium varies from 8.52 ppm in neat [C₄C₁im][PF₆] to 9.31 ppm in neat [C₄C₁im][MeSO₄] at room temperature,¹⁴ while for [C₈C₁im] cation the same resonance changes from 9.9 ppm in the chloride salt to 7.4 ppm in the tris(pentafluoroethyl)trifluorophosphate salt.¹⁵ Such significant environmental effects on the chemical shift are attributed to the difference in average bulk structure and chemical composition of the solvation shell of the reference cation. It is, therefore, still a structural effect on the chemical shift, but this time the “structural space” that one needs to explore in order to reproduce and eventually predict the chemical shift, is continuous rather than discrete. In fact, since it is necessary to account for the dynamics of the system, and since this is normally done using classical force fields (FFs) in a molecular dynamics (MD) simulation, a key ingredient to generate a reference dynamic trajectory is the choice of appropriate parameters for the interaction potential, including the internal bonding parameters, contact distance and well depth of the Lennard-Jones interaction, and partial charges of the atoms. To some extent, these parameters can be varied continuously, and thus resulting in a continuum spectrum of average bulk geometries of the ions in the bulk that can be generated by a classical MD simulation. It is noteworthy to mention that in such bulk dynamic systems, temperature effects can be expected to have an influence on the chemical shift as well. However, extensive studies on the temperature dependence of the NMR proton resonances of ILs are scarce. A notable exception is the recent work by Filippov and coworkers,¹⁶ where the Authors reported a non-negligible but rather weak temperature dependence for ethylammonium nitrate, with about 0.1 ppm of difference from 293 K up to 363 K. Analogous small effects are expected for ¹³C resonances, whose nuclei are less exposed to the environment and generally show a little dependence on the counter-anion.¹⁴ For these reasons, here we only focus on ¹H chemical shifts.

Some attempts have been made to predict the chemical shift of imidazolium ILs by DFT methods, and a comprehensive review of this subject can be found in ref. 13. It appears that the chemical shift of the cation's protons of imidazolium ILs with soft anions are much easier to predict than the analogous cases with hard anions, such as oxo-anions or chloride, where strong hydrogen bonds with the acidic ring protons, especially H2, occur in the bulk. In ref. 14 we reported the first computational NMR study of ILs, using a combined MD + DFT protocol, finding a nice agreement for imidazolium tetrafluoroborate salts (a soft anion) with MAE = 0.12 ppm, while for the methylsulfate salt (forming stronger hydrogen bonds through the negatively charged oxygen), the MAE was as large as 0.20 ppm. Similarly, some years ago, Chen and Izgorodina observed that imidazolium chloride and acetate salts also had larger average errors compared to other systems like tetrafluoroborate and hexafluorophosphate, when using a static cluster approach for the prediction of the chemical shift.¹⁷ More recently, concerning the difficulty of an accurate prediction of the chemical shift of neat [C₄C₁im]Cl, Lengvinaite et al. concluded that “Our results thus suggest that a refinement of the force field used in the present MD simulations of neat IL may be necessary in order to improve the local distribution of ions around the C2–H2 moiety of imidazolium cations”.¹⁸

In recent years, there has been a clear shift from using the so-called full-charge FFs,^19–24 where the sum of the atomic partial charges in a unit ion equals to ±1e, to the so-called scaled-charge FF²⁵ where the partial charges are proportionally reduced to account for the missing polarizability and charge transfer effect that takes place in the bulk phase of ILs.^26–28 Clearly, a more appropriate approach would be to use polarizable FFs,^29–32 which are, however, much more computationally demanding and, since now, they have not found a wide-spread use. Nonetheless, the research for more sophisticated FF for the simulation of ILs is still an active field with recent developments constantly presented in the literature, see for example the so-called virtual site FF (VSIL) presented by the group of Acevedo.^33,34 The fact that new FFs for ILs are still under development is a clear indication of the difficulties related with their parameterization, which is mostly the result of the scarcity of enthalpy-of-vaporization data as well as the strong intermolecular interactions responsible for the nano-structuring of ILs.

In this work, we present the results of a computational NMR study of two ILs, [C₄C₁im][BF₄] and [C₄C₁im]Cl, by using DFT calculations of the chemical shifts of the 1-butyl-3-methylimidazolium cation in clusters extracted from the trajectories obtained using four different FFs: a full-charge FF, a scaled-charge FF, a polarizable FF, and the recently derived VSIL FF mentioned above. We will show that the comparison of calculated and experimental data improves systematically as the quality of the FF is augmented, which for the first time enables one to obtain almost quantitative results concerning the chemical shift of imidazolium chloride ILs. The agreement of calculated and experimental results can also be used as a guide to support or discard microscopic structures obtained by various FFs.

Computational details

MD simulations

The structure, ring atoms numbering, and labeling of the alkyl groups of the 1-butyl-3-methylimidazolium cation as well as the tetrafluoroborate anion are shown in Chart 1.

Chart 1 Systems studied in this work. (a) Structural formula of the studied ions: 1-butyl-3-methylimidazolium, tetrafluoroborate, and chloride. For 1-butyl-3-methylimidazolium, we also show the ring atom numbering and labels (in italics) for the alkylic residues. (b) 1-butyl-3-methylimidazolium with the virtual-site charge in the VSIL FF shown as a dark dot placed in the plane of the imidazolium ring.

Four different FFs were tested for each ionic liquid: (1) a full-charge AMBER FF;³⁵ (2) a scaled-charge AMBER FF, where the partial charge on each atom is scaled by a factor of 0.807;³⁶ (3) a polarizable FF,^30,37 where the Thole's smearing dipole model with a factor of 2.392 is applied;^38,39 (4) an OPLS-VSIL FF,^33,34 where the partial charge on each ion is set to ±0.8e and an additional negative-partial-charge virtual site is placed inside each imidazolium ring, as shown in Chart 1(b). For the first two cases, the parameters are obtained from the model developed by Liu et al.⁴⁰ based on the AMBER FFs. For the polarizable FF, all the force parameters share the same values with the full-charge AMBER FF, except an extra polarizability for each atom.⁴¹ An iteration procedure is employed to solve the induced dipole on each atom.^37,42 Finally, the parameters of the OPLS-VSIL FF can be found in ref. 33 and 34.

All the relaxation and production simulations were performed with an isotropic barostat. In the simulations with the full-charge and scaled-charge AMBER FFs, the temperature and pressure were kept constant by using the Nosé–Hoover thermostat^43,44 with a time constant of 0.5 ps and the Parrinello–Rahman barostat^45,46 with a time constant of 1.0 ps, respectively. In the simulations with the polarizable FF, the barostat was changed to Nosé-Hoover^47–49 with the same time constant of 1.0 ps. The Nosé-Hoover thermostat and the Parrinello–Rahman barostat fail to control the temperature and pressure in the simulations of [C₄C₁im][BF₄] with the OPLS-VSIL FF, resulting in extremely large fluctuations for thermal quantities and a very slow diffusion. Therefore, the simulations with the OPLS-VSIL FF for both [C₄C₁im][BF₄] and [C₄C₁im]Cl alternatively employed the velocity rescaling thermostat⁵⁰ and the Berendsen barostat,⁵¹ both with a time constant of 1.0 ps, as reported in the original paper developing this FF.^33,34 The time integration was performed by the leap-frog algorithm with a timestep of 1 fs, except for the polarizable FF where a smaller timestep of 0.4 fs was used, as suggested in the previous study.³⁷ Simulations with the non-polarizable FFs were performed with Gromacs,⁵² and those with the polarizable FF were performed by a modified version of DL_POLY 2.^30,53,54

The simulated system contains 250 ion pairs in a cubic box, and periodic boundary conditions are applied on all three dimensions. The cutoff distances of both the van der Waals and the electrostatic interaction in the real space are set to be 13 Å. In the simulations with either the full-charge or the scaled-charge AMBER FF, starting from a random configuration equilibrated at T = 1000 K in the NVT ensemble for 2 ns, a simulated annealing procedure in the NPT ensemble (P = 1 atm) was performed thereafter with eight sequential steps: 2 ns at 1000 K, 2 ns at 800 K, 2 ns at 600 K, 2 ns at 550 K, 2 ns at 500 K, 2 ns at 450 K, 4 ns at 400 K, and 4 ns at 300 K. The final snapshot at T = 400 K with the full-charge FF was further equilibrated for 40 ns, serving as the initial configuration for both simulations with the polarizable FF and the OPLS-VSIL FF. Here we use the final configuration at T = 400 K since its dynamics is more flexible, allowing the system to effectively relax when changing the FF. In the simulations with the polarizable FF, the initial configuration was further equilibrated at T = 300 K for 20 ps before the production run, while in the simulations with the OPLS-VSIL FF, the initial configuration was relaxed at T = 300 K for 4 ns before sampling. The sample interval was set to be 10 ps for all the three non-polarizable FFs, resulting in 4000 configurations along a trajectory of 40 ns. In the later DFT calculations, only 400 of them with an even interval of 100 ps were used. The simulations with the polarizable FF are much slower, so a smaller sample interval of 2 ps was used, producing 400 configurations along a trajectory of 0.8 ns. The visualization of snapshots was realized with the VMD software,⁵⁵ while the analysis of trajectories was done with the GROMACS built-in routines, along with the TRAVIS software⁵⁶ to draw the spatial distribution functions.

DFT calculations

All DFT-NMR calculations were run at the GIAO⁵⁷ (Gauge-Independent Atomic Orbitals) B3LYP level⁵⁸ with the software package Gaussian 16.⁵⁹ This functional has been widely tested in the literature and it provides accurate chemical shifts for organic molecules.^1,60,61 Indeed, several other functionals (and methods, e.g. IGLO,⁶² Individual Gauge for Localized Orbitals; or CGST,⁶³ Continuous Set of Gauge Transformations) could have been considered, but the aim of this work is to evaluate trends as we change the average bulk geometry as obtained from different FFs, therefore just one reliable and standard protocol widely tested in the literature is sufficient. Several computational protocols have been used and they can be grouped into two sets. The first set consists of clusters extracted from the MD trajectory made by a reference [C₄C₁im] cation, selected at random for each configuration, plus 0, 1, 2, and 3 closest anions with respect to the geometric center of the ring. These ions are treated by DFT using the cc-pVTZ basis set.⁶⁴ Several basis sets for NMR have been presented in the literature, e.g. pcS-n⁶⁵ and pcSseg-n⁶⁶ by Jensen as well as pecS-n^67,68 by Rusakov and Rusakova. Nevertheless, the cc-pVTZ of triple-ζ quality basis set have been found to be very accurate for ¹H chemical shift calculations, and successfully used by ourselves in several previous studies.^60,61,69

The rest of the ions of the cluster are included in the DFT calculations just as point charges up to a cutoff radius of 20.0 Å, which is, in all cases, slightly less than half of the side length of the simulation box. Point charges are the same as used in the FFs described above (dipoles are not included in the polarizable FF). These protocols will be labeled hereafter as fcN, scN, polN, and vsilN, where N is the number of anions explicitly included (0, 1, 2, 3) and fc, sc, pol, and vsil indicate the full-charge, scaled-charge, polarizable and virtual-site FFs, respectively. The second set of protocols consists in extracting clusters with the same reference cation as above, but including a larger number of cations and anions. In particular we have included all ions which have at least one atom within 7.5 Å from the geometric center of the imidazolium ring. This choice ensures the inclusion of the entire first solvation shell, as evidenced by the radial distribution functions discussed below. In these cases, an ONIOM approach for the calculation of the NMR parameters is applied.⁷⁰ The reference cation is included in the High Level (HL), while the rest of ions are in the Low Level (LL). For one protocol, labeled hereafter as A, the HL basis set is cc-pVTZ while the LL basis set is 3-21g. For the other protocol noted as B, the HL and LL basis sets are 6-311g** and 6-31g*, respectively. Depending on the FF used to generate the trajectory, these protocols are labeled as fcA, scA, polA, and vsilA, along with fcB, scB, polB, and vsilB. The clusters for protocols A or B contained on average about 150 atoms, with only 25 atoms of the central [C₄C₁im] in the HL.

Finally, for all protocols, each proton-shielding constant is obtained as an average over the selected 400 configurations during the MD trajectory. A summary of all MD + DFT protocols is shown in Table 1 for reference. As an example, we show in Fig. 1 the most relevant type of clusters used.

Table 1 List of all the MD + DFT computational protocols used in this work. HL denotes the High Level ONIOM, while LL denotes the Low Level ONIOM

Label	MD FF	DFT level of theory	Cluster type
fc0	Full-charge	B3LYP/cc-pVTZ	Reference cation DFT, up to 20.0 Å point charges
fc1	Full-charge	B3LYP/cc-pVTZ	Reference cation + 1 anion DFT, up to 20.0 Å point charges
fc2	Full-charge	B3LYP/cc-pVTZ	Reference cation + 2 anion DFT, up to 20.0 Å point charges
fc3	Full-charge	B3LYP/cc-pVTZ	Reference cation + 3 anion DFT, up to 20.0 Å point charges
fcA	Full-charge	ONIOM HL = B3LYP/cc-pVTZ; LL = B3LYP/3-21g	Reference cation in HL, remaining ions with at least one atom within 7.5 Å in LL
fcB	Full-charge	ONIOM HL = B3LYP/6-311g**; LL = B3LYP/6-31g*	Reference cation in HL, remaining ions with at least one atom within 7.5 Å in LL
sc0	Scaled-charge	B3LYP/cc-pVTZ	Reference cation DFT, up to 20.0 Å point charges
sc1	Scaled-charge	B3LYP/cc-pVTZ	Reference cation + 1 anion DFT, up to 20.0 Å point charges
sc2	Scaled-charge	B3LYP/cc-pVTZ	Reference cation + 2 anion DFT, up to 20.0 Å point charges
sc3	Scaled-charge	B3LYP/cc-pVTZ	Reference cation + 3 anion DFT, up to 20.0 Å point charges
scA	Scaled-charge	ONIOM HL = B3LYP/cc-pVTZ; LL = B3LYP/3-21g	Reference cation in HL, remaining ions with at least one atom within 7.5 Å in LL
scB	Scaled-charge	ONIOM HL = B3LYP/6-311g**; LL = B3LYP/6-31g*	Reference cation in HL, remaining ions with at least one atom within 7.5 Å in LL
pol3	Polarizable	B3LYP/cc-pVTZ	Reference cation + 3 anion DFT, up to 20.0 Å point charges
polA	Polarizable	ONIOM HL = B3LYP/cc-pVTZ; LL = B3LYP/3-21g	Reference cation in HL, remaining ions with at least one atom within 7.5 Å in LL
vsil3	Virtual-site	B3LYP/cc-pVTZ	Reference cation + 3 anion DFT, up to 20.0 Å point charges
vsilA	Virtual-site	ONIOM HL = B3LYP/cc-pVTZ; LL = B3LYP/3-21g	Reference cation in HL, remaining ions with at least one atom within 7.5 Å in LL

---

Fig. 1 Example of clusters extracted from the MD trajectory with the full-charge FF. All clusters include the reference [C₄C₁im] cation at the center, for which the proton chemical shifts are calculated. In addition, (a) and (b) include the three closest anions of BF₄⁻ (a) or Cl⁻ (b), to be treated at the DFT level, and the rest of the ions up to 20.0 Å from the geometric center of the imidazolium ring as point charges (represented as red dots). Protocols fc3, sc3, pol3, and vsil3 are applied to these types of clusters. (c) and (d) include all the ions (either positive or negative) for which at least one atom is found within 7.5 Å from the geometric center of the imidazolium ring. The central reference cation is treated at the ONIOM high level, while the rest of the ions are treated at the ONIOM low level. Protocols fcA/B, scA/B, polaA, and vsilA are applied to these clusters separately.

Not all protocols have been applied to all systems: for testing purposes, we applied protocols fc0, fc1, fc2, fc3, fcA, fcB, and sc0, sc1, sc2, sc3, scA, scB to both salts. By contrast, for the polarizable FF and VSIL FF, we only applied the “best” protocols polA, vsilA, pol3, and vsil3. The total number of DFT calculations amounts to 12 800 (16 protocols × 2 salts × 400 configurations).

Experimental ¹H chemical shifts have been taken from ref. 18 and 71 for the chloride salt and ref. 72 for [BF₄] salt, which are reported in Table S1 in SI. In order to compare the calculated and experimental chemical shifts, we first rescaled all data with respect to the resonance (either experimental or calculated) of the terminal methyl group of the butyl chain, which is the most shielded and least affected one. This procedure is very useful to avoid the calculation of the shielding constant of the NMR reference (tetramethylsilane) using the same protocol, which would otherwise require a MD simulation of TMS in a given solvent modeled with an appropriate FF.^14,18 Moreover, it provides more accurate results, and thus enhances the predictive power of the computational protocol.⁷³ Therefore, the rescaled chemical shifts are calculated as

δ_expt,resc = δ_expt − δ_expt(CH₃−1)

and

δ_calc,resc = σ_calc(CH₃−1) − σ_calc,

where δ is the chemical shift and σ the shielding constant. Unless otherwise stated, δ_calc and δ_expt will refer to the rescaled parameters from now on. We then calculated the slope a, intercept b, and correlation coefficient R² of the linear fitting,

δ_calc = aδ_expt + b,

the maximum absolute error, MaxErr, defined as

MaxErr = max(|δ_calc − δ_expt|)

and the mean absolute error, MAE, defined as

In all these statistical evaluations, N represents the number of resonances in the ¹H spectrum except the data point of the terminal methyl group, which is used as a reference and therefore excluded. These parameters are finally used to judge the degree of agreement between calculated and experimental data.

Results and discussion

MD simulations

We first take a look at the spatial correlation between different ions reflected by the center-of-mass partial radial distribution functions (PRDFs) of cation–cation g₊₊(r), anion–anion g₋₋(r), and cation–anion g₊₋(r).

As shown in Fig. 2, for both [C₄C₁im]Cl and [C₄C₁im][BF₄], the full-charge AMBER, scaled-charge AMBER, and polarizable FFs produce almost the same structures on large scales. However, a qualitatively different local structure for the second shell of the anion, corresponding to the first peak of g₋₋(r) in Fig. 2(c) and (f), is exhibited in the simulation results with the polarizable and non-polarizable AMBER FFs, consistent with the previous study on another kind of ionic liquid.^30,37 The second shell of the cation, corresponding to the first peak of the g₊₊(r) in Fig. 2(a) and (d), also exhibits a difference in distance that the polarizable FF indicates a slightly closer one: In the system of [C₄C₁im]Cl, 0.81 nm for the full-charge AMBER FF, 0.83 nm for the scaled-charge FF, while 0.735 nm for the polarizable FF; in the system of [C₄C₁im][BF₄], 0.8 nm for both full-charge and scaled-charge AMBER FFs, while 0.755 nm for the polarizable FF. The correlation between cation and anion produced by these three different FFs are quite similar, evidenced by the nearly overlapping curves of g₊₋(r) shown in Fig. 2(b) and (e). The simulations with the OPLS-VSIL FF produce qualitatively different results for both g₊₋(r) and g₊₊(r), demonstrating an obviously closer distance of both the first shell (anion) and the second shell (cation) around a cation. In the system of [C₄C₁im]Cl, the distances are 0.41 nm and 0.69 nm, respectively, smaller than 0.445 nm and 0.735 nm from the polarizable FF; in the system of [C₄C₁im][BF₄], the ones from OPLS-VSIL are 0.45 nm and 0.73 nm, also smaller than the ones from the polarizable FF, 0.455 nm and 0.755 nm. It should be noticed that this tendency can be extended to a distance as far as 1.5 nm. These differences are more obvious for the case of [C₄C₁im]Cl than [C₄C₁im][BF₄]. The first peak of g₋₋(r) for [C₄C₁im][BF₄] gained from OPLS-VSIL, as depicted in Fig. 2(f), indicates little deviation from the AMBER FF. However, as shown in Fig. 2(c), OPLS-VSIL FF gives closer peaks than the other three FFs, except the first one in g₋₋(r) for [C₄C₁im]Cl. The stronger attraction created by the OPLS-VSIL FF is also evidenced by a higher value of the first peak of all the three PRDFs, manifesting a more ordered local structure.

Fig. 2 Partial radial distribution functions of (a)–(c) [C₄C₁im]Cl and (d)–(f) [C₄C₁im][BF₄]. (a) and (d) Cation–cation. (b) and (e) Cation–anion. (c) and (f) Anion–anion. It can be seen that different force fields can affect the local structure characterized by the first peak, while only the influence of OPLS-VSIL can extend to a relatively large scale. Moreover, the differences in [C₄C₁im]Cl are generally larger than the ones in [C₄C₁im][BF₄].

A detailed description of local structures can be obtained by the spatial distribution function (SDF). In Fig. 3, we show the SDF of anions around a central cation in [C₄C₁im]Cl. It can be seen that the iso-surface has a relatively large hole along the direction of the C2–H2 bond for the full-charge AMBER, the scaled-charge AMBER, and the polarizable FF, indicating a local minimum of probability and thus a relatively high energy barrier there. However, this hole does not appear with the OPLS-VSIL FF, in agreement with the previous AIMD results,³⁴ as well as with the results of Car–Parrinello MD simulations of the analogous system [C₂C₁im]Cl.^74,75 On the other hand, as shown in Fig. 4, the SDFs of [BF₄] around [C₄C₁im] obtained with the first three FFs are also qualitatively the same, while a clear difference, similar to the case of the chloride salt, is observed for the VSIL FF.

Fig. 3 Spatial density function of the probability to find the chloride anion at a given position around the imidazolium ring. The iso-density value is set to be 20. (a) Full-charge FF. (b) Scaled charge FF. (c) Polarizable FF. (d) VSIL FF.

Fig. 4 Spatial density function of the probability to find the boron atom of [BF₄] at a given position around the imidazolium ring. The iso-density value is set to be 20. (a) Full-charge FF. (b) Scaled-charge FF. (c) Polarizable FF. (d) VSIL FF.

It is worth mentioning that the SDF provides a qualitative rather than quantitative characterization of the local structure, where the shape of the iso-surface highly depends on the iso-value. For instance, we can see a hole on the iso-surface along the C2–H2 direction in the [C₄C₁im]Cl system with the VSIL FF when the iso-value is larger than 30, manifesting the existence of a local minimum of probability connecting the two high-probability lobes in the SDF.

DFT calculations

The statistical parameters of the correlation between calculated and experimental chemical shifts are summarized in Tables 2 and 3. Additional figures are reported in SI.

Table 2 Correlation parameters for various computational protocols applied to the [C₄C₁im][BF₄] system. MAE refers to mean absolute error, MaxErr to maximum error, Err H2 to absolute error for proton H2, R² to correlation coefficient of the linear fit, a to the slope of the linear fitting line, and b to the intercept of the linear fitting line. Uncertainties are: 0.005 ppm on the MAE; 0.05 ppm on the intercept; 0.01 on the slope

	MAE/ppm	MaxErr/ppm	Err H2/ppm	R ^2	a	b/ppm
fcA	0.25	0.52	0.14	0.9989	1.0272	0.13
fcB	0.16	0.38	0.09	0.9985	1.0106	0.07
fc0	0.17	0.56	0.56	0.9983	0.9673	−0.02
fc1	0.14	0.42	0.42	0.9989	0.9781	−0.03
fc2	0.11	0.29	0.29	0.9992	0.9916	−0.04
fc3	0.09	0.21	0.21	0.9991	1.0027	−0.05
scA	0.19	0.48	0.01	0.9986	1.0208	0.10
scB	0.14	0.35	0.20	0.9982	1.0051	0.05
sc0	0.19	0.64	0.64	0.9979	0.9700	−0.06
sc1	0.15	0.50	0.50	0.9987	0.9786	−0.06
sc2	0.11	0.37	0.37	0.9990	0.9912	−0.06
sc3	0.11	0.26	0.26	0.9990	1.0045	−0.06
polA	0.22	0.44	0.21	0.9993	1.0385	0.05
pol3	0.08	0.25	0.02	0.9994	1.0171	−0.05
vsilA	0.19	0.39	0.17	0.9982	0.9865	0.19
vsil3	0.14	0.40	0.40	0.9992	0.9560	0.11

---

Table 3 Correlation parameters for various computational protocols applied to the [C₄C₁im]Cl system. MAE refers to mean absolute error, MaxErr to maximum error, Err H2 to absolute error for proton H2, R² to correlation coefficient of the linear fit, a to the slope of the linear fitting line, and b to the intercept of the linear fitting line. Uncertainties are: 0.01 ppm on the MAE; 0.1 ppm on the intercept; 0.02 on the slope

	MAE/ppm	MaxErr/ppm	Err H2/ppm	R ^2	a	b/ppm
fcA	0.34	1.25	1.25	0.9961	0.8917	0.24
fcB	0.39	1.36	1.36	0.9959	0.8871	0.20
fc0	0.80	2.23	2.23	0.9950	0.8145	0.12
fc1	0.68	1.79	1.79	0.9978	0.8452	0.08
fc2	0.56	1.47	1.47	0.9983	0.8769	0.04
fc3	0.45	1.28	1.28	0.9980	0.9016	0.03
scA	0.51	1.53	1.53	0.9972	0.8622	0.21
scB	0.63	1.81	1.81	0.9957	0.8496	0.14
sc0	0.96	2.55	2.55	0.9941	0.7877	0.08
sc1	0.86	2.21	2.21	0.9964	0.8131	0.06
sc2	0.73	1.87	1.87	0.9976	0.8436	0.04
sc3	0.62	1.62	1.62	0.9977	0.8681	0.03
polA	0.44	1.14	1.14	0.9971	0.9098	0.05
pol3	0.46	0.99	0.99	0.9991	0.9149	−0.02
vsilA	0.31	0.91	0.91	0.9943	0.9494	0.09
vsil3	0.22	0.63	0.63	0.9976	0.9657	−0.01

---

First, we briefly discuss the case of tetrafluoroborate. The time evolution of the shielding constants for various protocols and imidazolium resonances can be found in SI, Fig. S3–S46. We do not discuss in detail these data, but it is noteworthy that the instantaneous fluctuations of the shielding constant are very large, so several configurations need to be considered to achieve a reasonable averaged value. Just as an example, as shown in Fig. S6 in SI, the shielding constant calculated with the fc3 protocol spans a range from a minimum of 21.12 ppm (absolute value without rescaling) to a maximum of 25.15 ppm, with an error of 0.036 ppm on the mean value after averaged over 400 configurations. These figures can be considered as representative of all protocols.

In Fig. 5(a), we show the correlation between experimental and calculated chemical shifts data (rescaled with respect to the resonance of the methyl group of the butyl chain) with protocols FF3, where FF is either fc, sc, pol or vsil. As we can see, all FFs perform quite well. Although there are small variations among different methods, a close inspection of Fig. 6(a) and (b) reveals that, as might be expected, increasing the number of anions close to the imidazolium ring explicitly included in the DFT calculation improves the overall agreement.

Fig. 5 Correlation between calculated and experimental chemical shifts of [C₄C₁im] for some selected computational protocols in the tetrafluoroborate salt (a) and the chloride salt (b). The values are rescaled with respect to the shielding constant of the methyl group of the butyl chain.

Fig. 6 Statistical parameters of the correlation between experimental and calculated [C₄C₁im] chemical shifts. (a) Mean absolute error in [C₄C₁im][BF₄]. (b) Maximum error in [C₄C₁im][BF₄]. (c) Mean absolute error in [C₄C₁im]Cl. (d) Maximum error in [C₄C₁im]Cl.

With the pol3 protocol, a MAE of just 0.08 ppm is obtained. Somewhat unexpectedly, the vsil3 protocol performance is slightly worse, although still exhibits relatively high accuracy. Besides, both the fc3 and sc3 protocols show very good performance. In short, the strategy of including the three closest anions to the ring in the DFT calculation and treating the rest of the bulk phase with point charges appears to work very well for predicting the chemical shift. Surprisingly, the other two types of protocols, A and B, based on the ONIOM scheme and using larger clusters at the DFT level, have somewhat worse performance irrespective of the FF used to generate the trajectories. However, as a general conclusion, we can say that a very accurate prediction of proton chemical-shifts can be achieved with all types of FFs using several computational approaches, in agreement with previous results already discussed in the Introduction section. This is an indication, on the one hand, of the correctness of the approach, but on the other hand, of a relatively poor sensitivity of the chemical shift of imidazolium cations on the environment in the case of soft anions like tetrafluoroborate. These results are consistent with the outcomes of a very recent study concerning imidazolium tetrafluoroborate and its mixtures with water.⁷⁶

The situation changes completely in the case of the chloride system. In Fig. 5(b) we show the analogous correlation between calculated and experimental chemical shifts already seen for the tetrafluoroborate salt. Clearly, the ring protons, especially H2, are significantly more in error when compared with the same protocol applied to the [BF₄] case, which can be seen from MAE and MaxErr for the chloride case shown in Fig. 6(c) and (d), respectively. First, all the values are generally much larger, with discrepancies in the MaxErr (which almost in all cases is due to H2) even larger than 2.5 ppm. The effect of including additional anions in the cluster is quite evident. For example, the MAE varies from 0.80 ppm to 0.45 ppm when the protocol goes from fc0 to fc3. Somewhat worse is the performance of the scaled-charge FF, compared to the full-charge FF, since the MAE changes from 0.96 to 0.62 ppm when the protocol goes from sc0 to sc3.

A comment is necessary concerning protocols A and B in the case of chloride salts (while no issues occurred with protocols 0–3). For protocols A and B, we noted a significant slower convergence rate of the self-consistent field (SCF) procedure during the electronic structure calculations. In a limited number of cases, the SCF had to be converged using the quadrature convergence method as implemented in Gaussian 16. These issues are discussed in details in SI, see Fig. S2.

At any rate, it appears more convenient to focus on protocols 0–3 with various FFs since no SCF convergence issues were detected. We have discussed already the poor results for the full-charge and scaled-charge FFs. Instead, for the polarizable FF, we note a non-negligible improvement: Although the MAE at the pol3 level is practically the same as for the fc3 level (0.46 ppm vs. 0.45 ppm), the MaxErr, which coincides with the error on proton H2, is 1.28 ppm and 1.62 ppm for fc3 and sc3, respectively, but it drops down to 1.14 ppm at the pol3 level. Finally, a significant improvement is obtained with the VSIL trajectory: here the MAE is just 0.22 ppm and the MaxErr 0.63 ppm, providing a very good result which brings the statistical performance of the computational protocol based on the MD + DFT vsil3 quite close to the performance obtained for the softer anions like tetrafluoroborate.

As can be judged from Fig. 5, various FFs in the case of the chloride salts fail to predict the extremely large deshielding of the ring protons, particularly H2. While the experimental value for H2 is a remarkable 10.23 ppm (rescaled value of 9.73 ppm w.r.t. the methyl group on the butyl chain), the computed values always lie below (rescaled calculated values of 8.45, 8.11, 8.74, and 9.10 ppm for fc3, sc3, pol3, and vsil3 protocols). This lack of success should be attributed to the poor ability of classical FFs in describing correctly the hydrogen bonds (HBs) formed by the proton H2 on the imidazolium ring and the chloride ion.

Potential energy surface

To gain further insight on the origin of the inaccuracy described above, we have calculated the potential energy surface (PES) of a model ion pair, especially concerning the position of the chloride with respect to the proton H2 of the imidazolium ring. The simpler and symmetric cation dimethylimidazolium [C₁C₁im] is used for this numerical calculation to avoid the diversity of conformations caused by the long side chain in [C₄C₁im]. The geometry of [C₁C₁im] is optimized at the ωB97XD/6-311+g** level.⁷⁷ The PES is then scanned in the 2D plane of the imidazolium ring as a function of the C2–Cl distance and the angle θ between the C2–H2 direction and the C2–Cl direction, see Fig. 7.

Fig. 7 Potential energy surface experienced by the chloride anion interacting with the dimethylimidazolium cation. (a) Schematic illustration of the molecular structure of dimethylimidazolium and the definition of r and θ. Here θ = 0 corresponds to the direction along the C2–H2 bond, with an anti-clockwise rotation corresponding to a positive value. (b) Full-charge AMBER force field. (c) Scaled-charge AMBER force field. (d) Polarizable force field. (e) OPLS-VSIL force field. (f) DFT calculation in solvent environment, see text for details. The color in each panel of (b)–(f) corresponds to the relative value of the potential energy, where the lowest energy of each case has been shifted to zero. The color bar applies at a logarithmic scale.

The calculations have been run at two different levels of theory: one using the four classical FFs for the MD trajectories and the other at the DFT level (ωB97XD(PCM)/6-311+g**) as a reference for the PES. The choice of the DFT level of theory, which includes the solvent effects (acetonitrile) by means of the Polarizable Continuum Model⁷⁸ to dump the electrostatic interactions in gas phase, guarantees that the long-range solvent effects and dispersive interactions are properly accounted for. These effects, empirically included in the FF by means of LJ interactions, play a prominent role in determining the structure of non-covalent pairs. On the other hand, they have a negligible impact on NMR properties.^14,17,74,79 For the classical FF protocols, the partial charges of the atoms on the side chains are copied from the same atomic group on [C₄C₁im], and the partial charges of all the atoms on the cation are so scaled that their total sum is 1 for the full-charge AMBER and polarizable FFs, 0.807 for the scaled-charge AMBER FF, and 0.8 for the OPLS-VSIL FF. The potential energy is calculated for the ion pair without including the intramolecular contributions of the cation and taking the van der Waals and electrostatic interactions into account for the non-polarizable FFs, and also the charge-induced dipole, induced dipole-induced dipole, and self-energy of the induced dipole for the polarizable FF. Although the polarization is a highly collective effect, here we only consider this effect for one ion pair, which possibly reduces the difference with the non-polarizable FFs. Due to the mirror-symmetry of the molecule and the planar geometry of the imidazolium ring, the PES is symmetric with respect to θ = 0. Here we focus on two key quantities reflected by the PES: the distance between Cl and H2 at the energy minima and the relative height of the energy barrier at θ = 0. The DFT result, which accurately accounts for hydrogen bonds as well as polarizability and possible charge-transfer effects, suggests a much flatter energy landscape, and especially a much lower energy barrier between the two symmetric energy basins compared to the one from empirical FFs. The potential energy landscapes with the full-charge AMBER and polarizable FFs are quite similar, the one with the scaled-charge AMBER FF slightly shifts to an even farther distance, and the one with the OPLS-VSIL FF shifts to a closer distance, but they all retain similar morphologies. Detailed data of the distance of H2⋯Cl and the height of energy barrier listed in Table 4 demonstrate that the full-charge FF predicts a distance for the minimum energy configuration shorter than the scaled-charge FF but slightly larger than the polarizable one, which are all significantly larger than the distance given by DFT calculation. The OPLS-VSIL FF, instead, perfectly predicts the distance but still shows a higher energy barrier compared to the DFT results.

Table 4 Energy barrier, ΔE, and the distance between the Cl and H2 atoms of the potential energy landscapes in Fig. 7 using various computational protocols

	Full-charge FF	Scaled-charge FF	Polarizable FF	OPLS-VSIL FF	DFT
ΔE/kJ mol^−1	6.60	3.81	7.28	3.93	0.91
R(Cl⋯H2)/Å	2.66	2.74	2.62	2.31	2.33

---

The above PESs suggest that the error stems from the fact that the closest distances between chloride and H2 given by the AMBER-based FFs are too large. Considering the fact that the key parameter is the Lennard-Jones contact distance σ of the ring hydrogen, the AMBER-based FFs all employ values typical of aromatic hydrogens, in our case 2.42146 Å for the full-charge, scaled-charge and polarizable FF, and 1.85 Å for the VSIL FF. Clearly, the smaller value for the VSIL FF is the main reason of the large improvement obtained using these trajectories. However, it is well-known that, to account for true HB in classical FF, e.g., with hydroxyl groups, a contact distance of zero is set for the hydrogen atoms. This workaround allows to account for the short O⋯O distances in hydrogen-bonded alcohols and/or water molecules in classical MD simulations. In fact, such strong HBs have a non-negligible covalent character, and they can only be properly described by quantum chemistry calculation.⁸⁰ The hydrogen bond formed between the ring protons and hard anions, such as chlorides, is certainly weaker than typical hydroxyl hydrogen bonds. Therefore, setting the contact distance to zero for the ring protons of the imidazolium ring would not be a proper solution. We should recall that this issue affects almost all of the available FFs for ILs presented in the literature. For example, the popular CL&P family developed by Canongia-Lopes and Padua also features an LJ contact distance parameter σ of the ring protons of 2.42 Å,²¹ practically the same value as in the AMBER FF used here. Therefore, since this is the most important parameter influencing the H-bond, we expect similar results to be obtained.

A similar qualitative result was obtained several years ago by comparing the predictive power of a classical full-charge FF with that of a Car–Parrinello simulation (CPMD) for ethylmethylimidazolium chloride.⁷⁴ However, the limited sampling in time and size of the CPMD run available at that time did not allowed an thorough comparison. In the PESs reported in Fig. 7, all FFs clearly exhibit two symmetric energy minima with a relatively high barrier to pass from one basin to the other, which corresponds to a minimum of probability to find the chloride anion just in front of H2. By contrast, the DFT PES does show an almost barrierless profile connecting the two minima, in agreement with the SDF reported in the literature.^74,75 We believe, therefore, that two key elements should be taken into account when developing the next type of classical FF for accurate simulations of imidazolium chlorides. On the one hand, the set of parameters should correctly describe the distance of minimum approach between chloride and imidazolium ring protons, and on the other hand, they need to provide an almost flat PES in front of the H2 ring proton.

Conclusions

In summary, we have investigated the computational NMR of two kinds of ionic liquids: [C₄C₁im][BF₄] and [C₄C₁im]Cl, providing a detailed comparison among different classical FFs as well as different DFT protocols. The results highlight the significant effects of the hydrogen bond, which is usually poorly described by classical FFs: for a softer anion BF₄ with weaker hydrogen bond, all the classical FFs perfectly predict the chemical shift value consistent with the experiments; by contrast, for Cl with stronger hydrogen bond effect, a significant improvement is observed from full-charge FF to polarizable FF to VSIL one. A mechanistic analysis by considering the PES produced by different FFs reveals two essential elements for an accurate local structure: the position of the local minimum and the barrier between two symmetric minima. It should be noticed that even VSIL predicts a significantly higher energy barrier, leaving some space for future modification.

Due to the limited ability of experiments to detect the structures on single molecular scale, especially dynamically and in complex environment, along with the limited system size able to be reached by the ab-initio calculation, molecular dynamics simulation plays a critical role in the study of ionic liquids. Consequently, a high-quality force field is indispensable and fundamental for this area. Previous researchers made great contributions to this topic from several aspects, e.g., the scaled-charge FF includes a scale factor on the partial charge aimed at a better description on the diffusion constants, while the polarizable and the VSIL FFs focus on the local structures compared with the ones from ab initio calculation. Computational NMR is known to be sensitive to the local environment, and therefore can be used as a fingerprint of the local structures. With this powerful tool, we are able to offer a quantitative rather than qualitative criterion to tell which force field gives a more realistic local structure as well as point out how far we are still away from “accurate”. This work is anticipated to not only provide insights for the verification of the classical FFs, but also be heuristic for the development of the next generation FFs.

Author contributions

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: Shielding constants as a function of time for all protocols used. See DOI: https://doi.org/10.1039/d5cp03322e.

Acknowledgements

We thank Mr Tao Wang for his kind help on the calculation of spatial distribution functions. DFT calculations were run using several computational facilities: the cluster of the C3P community in Padova (https://www.chimica.unipd.it/en/node/319), the cluster of the CloudVeneto infrastructure in Padova (https://cloudveneto.it/), the HPC of CINECA Consortium in Bologna (https://www.cineca.it) thanks to ISCRA projects HP10CQB2Q7 and HP10CM7DEJ, and the HPC cluster of ITP-CAS. We thank all these institutions for the allocation of computational time. This work is financially supported by the CNR-CAS 2024-2026 bilateral agreement titled “Membranes for water desalination based on discotic ionic liquid crystals”.

Notes and references

1. A. Bagno, Chem. – Eur. J., 2001, 7, 1652–1661.
2. G. Barone, L. Gomez-Paloma, D. Duca, A. Silvestri, R. Riccio and G. Bifulco, Chem. – Eur. J., 2002, 8, 3233–3239.
3. J. Meiler, E. Sanli, J. Junker, R. Meusinger, T. Lindel, M. Will, W. Maier and M. Köck, J. Chem. Inf. Comput. Sci., 2002, 42, 241–248.
4. D. A. Forsyth and A. B. Sebag, J. Am. Chem. Soc., 1997, 119, 9483–9494.
5. S. G. Smith and J. M. Goodman, J. Am. Chem. Soc., 2010, 132, 12946–12959.
6. S. G. Smith and J. M. Goodman, J. Org. Chem., 2009, 74, 4597–4607.
7. M. O. Marcarino, S. Cicetti, M. M. Zanardi and A. M. Sarotti, Nat. Prod. Rep., 2022, 39, 58–76.
8. M. O. Marcarino, M. M. Zanardi, S. Cicetti and A. M. Sarotti, Acc. Chem. Res., 2020, 53, 1922–1932.
9. A. M. Sarotti, Org. Biomol. Chem., 2013, 11, 4847–4859.
10. F. L. P. Costa, A. C. F. de Albuquerque, R. G. Fiorot, L. M. Lião, L. H. Martorano, G. V. S. Mota, A. L. Valverde, J. W. M. Carneiro and F. M. dos Santos Junior, Org. Chem. Front., 2021, 8, 2019–2058.
11. A. C. F. de Albuquerque, L. H. Martorano and F. M. dos Santos, Front. Nat. Prod., 2024, 2, 1321043.
12. G. Saielli, Adv. Theory Simulations, 2018, 1, 1800084.
13. G. Saielli, NMR Spectroscopic Parameters: Theories and Models, Computational Codes and Calculations, ed G. A. Aucar, RSC Publishing, London, 2025, pp. 395–429.
14. A. Bagno, F. D’Amico and G. Saielli, J. Phys. Chem. B, 2006, 110, 23004–23006.
15. T. Cremer, C. Kolbeck, K. R. J. Lovelock, N. Paape, R. Wölfel, P. S. Schulz, P. Wasserscheid, H. Weber, J. Thar, B. Kirchner, F. Maier and H.-P. Steinrück, Chem. – Eur. J., 2010, 16, 9018–9033.
16. O. I. Gnezdilov, O. N. Antzutkin, R. Gimatdinov and A. Filippov, Magn. Reson. Imaging, 2020, 74, 84–89.
17. S. Chen and E. I. Izgorodina, Phys. Chem. Chem. Phys., 2017, 19, 17411–17425.
18. D. Lengvinaitė, S. Kvedaraviciute, S. Bielskutė, V. Klimavicius, V. Balevicius, F. Mocci, A. Laaksonen and K. Aidas, J. Phys. Chem. B, 2021, 125, 13255–13266.
19. Y. Wang and G. A. Voth, J. Am. Chem. Soc., 2005, 127, 12192–12193.
20. J. C. Lopes, J. Deschamps and A. A. H. Padua, J. Phys. Chem. B, 2004, 108, 11250.
21. J. N. Canongia Lopes, J. Deschamps and A. A. H. Pádua, J. Phys. Chem. B, 2004, 108, 2038–2047.
22. J. N. A. Canongia Lopes and A. A. H. Pàdua, J. Phys. Chem. B, 2006, 110, 3330–3335.
23. J. N. Canongia Lopes and A. A. H. Pádua, Theor. Chem. Acc., 2012, 131, 1129.
24. T. Köddermann, D. Paschek and R. Ludwig, Chem. Phys. Chem., 2007, 8, 2464–2470.
25. Z. Sun, Z. Gong, L. Zheng, P. Kalhor, Z. Huai and Z. Liu, J. Ion. Liq., 2022, 2, 100043.
26. F. Philippi, K. Goloviznina, Z. Gong, S. Gehrke, B. Kirchner, A. A. H. Pádua and P. A. Hunt, Phys. Chem. Chem. Phys., 2022, 24, 3144–3162.
27. C. Schroder, Phys. Chem. Chem. Phys., 2012, 14, 3089–3102.
28. Y. Zhang and E. J. Maginn, J. Phys. Chem. Lett., 2015, 6, 700–705.
29. K. Goloviznina, J. N. Canongia Lopes, M. Costa Gomes and A. A. H. Pádua, J. Chem. Theory Comput., 2019, 15, 5858–5871.
30. T. Yan, C. J. Burnham, M. G. Del Pòpolo and G. A. Voth, J. Phys. Chem. B, 2004, 108, 11877–11881.
31. C. Schröder, A. Lyons and S. W. Rick, Phys. Chem. Chem. Phys., 2020, 22, 467–477.
32. N. Wang and E. J. Maginn, J. Phys. Chem. B, 2024, 128, 871–881.
33. B. Doherty, X. Zhong and O. Acevedo, J. Phys. Chem. B, 2018, 122, 2962–2974.
34. K. Yue, B. Doherty and O. Acevedo, J. Phys. Chem. B, 2022, 126, 3908–3919.
35. W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell and P. A. Kollman, J. Am. Chem. Soc., 1995, 117, 5179–5197.
36. V. V. Chaban, I. V. Voroshylova and O. N. Kalugin, Phys. Chem. Chem. Phys., 2011, 13, 7910–7920.
37. T. Yan, Y. Wang and C. Knox, J. Phys. Chem. B, 2010, 114, 6905–6921.
38. B. T. Thole, Chem. Phys., 1981, 59, 341–350.
39. C. J. Burnham, J. Li, S. S. Xantheas and M. Leslie, J. Chem. Phys., 1999, 110, 4566–4581.
40. Z. Liu, S. Huang and W. Wang, J. Phys. Chem. B, 2004, 108, 12978–12989.
41. Y. Gu and T. Yan, J. Phys. Chem. A, 2013, 117, 219–227.
42. J. Applequist, J. R. Carl and K.-K. Fung, J. Am. Chem. Soc., 1972, 94, 2952–2960.
43. S. Nosé, Mol. Phys., 1984, 52, 255–268.
44. W. G. Hoover, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 31, 1695–1697.
45. M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182–7190.
46. S. Nosé and M. L. Klein, Mol. Phys., 1983, 50, 1055–1076.
47. W. G. Hoover, Phys. Rev. A: At., Mol., Opt. Phys., 1986, 34, 2499–2500.
48. G. J. Martyna, D. J. Tobias and M. L. Klein, J. Chem. Phys., 1994, 101, 4177–4189.
49. W. Shinoda, M. Shiga and M. Mikami, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 134103.
50. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 14101.
51. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690.
52. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1–2, 19–25.
53. W. Smith, Comput. Phys. Commun., 1992, 67, 392–406.
54. W. Smith and T. R. Forester, J. Mol. Graph., 1996, 14, 136–141.
55. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graph., 1996, 14, 33–38.
56. M. Brehm, M. Thomas, S. Gehrke and B. Kirchner, J. Chem. Phys., 2020, 152, 164105.
57. K. Wolinski, J. F. Hinton and P. Pulay, J. Am. Chem. Soc., 1990, 112, 8251–8260.
58. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
59. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. J. Montgomery, J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16, Revision B.01, Gaussian, Inc., Wallingford CT, 2016.
60. A. Bagno, F. Rastrelli and G. Saielli, J. Phys. Chem. A, 2003, 107, 9964–9973.
61. A. Bagno, F. Rastrelli and G. Saielli, Chem. – Eur. J., 2006, 12, 5514–5525.
62. W. Kutzelnigg, U. Fleischer and M. Schindler, The IGLO-Method: Ab-initio Calculation and Interpretation of NMR Chemical Shifts and Magnetic Susceptibilities, ed. U. Fleischer, W. Kutzelnigg, H.-H. Limbach, G. J. Martin, M. L. Martin and M. Schindler, Springer Berlin Heidelberg, Berlin, Heidelberg, 1991, pp. 165–262.
63. T. A. Keith and R. F. W. Bader, Chem. Phys. Lett., 1993, 210, 223–231.
64. T. H. Dunning, J. Chem. Phys., 1989, 90, 1007–1023.
65. F. Jensen, J. Chem. Theory Comput., 2008, 4, 719–727.
66. F. Jensen, J. Chem. Theory Comput., 2015, 11, 132–138.
67. Y. Y. Rusakov and I. L. Rusakova, Phys. Chem. Chem. Phys., 2023, 25, 18728–18741.
68. Y. Y. Rusakov and I. L. Rusakova, J. Chem. Phys., 2022, 156, 244112.
69. S. G. G. G. M. I. B. A. Tähtinen, Chem. – Eur. J., 2008, 14, 10445–10452.
70. P. B. Karadakov and K. Morokuma, Chem. Phys. Lett., 2000, 317, 589–596.
71. S. Cha, M. Ao, W. Sung, B. Moon, B. Ahlström, P. Johansson, Y. Ouchi and D. Kim, Phys. Chem. Chem. Phys., 2014, 16, 9591–9601.
72. P. A. Z. Suarez, S. Einloft, J. E. L. Dullius, R. F. de Souza and J. Dupont, J. Chim. Phys., 1998, 95, 1626–1639.
73. J. Wang, Y. Liu, W. Li and G. Gao, RSC Adv., 2018, 8, 28604–28612.
74. A. Bagno, F. D’Amico and G. Saielli, ChemPhysChem, 2007, 8, 873–881.
75. B. L. Bhargava and S. Balasubramanian, Chem. Phys. Lett., 2006, 417, 486–491.
76. E. Sipavičius, G. Majauskaitė, D. Lengvinaitė, S. Bielskutė, V. Klimavicius, F. Mocci, A. Laaksonen and K. Aidas, Chem. Phys. Lett., 2025, 877, 142295.
77. J.-D. Chai and M. Head-Gordon, Phys. Chem. Chem. Phys., 2008, 10, 6615–6620.
78. J. Tomasi, B. Mennucci and R. Cammi, Chem. Rev., 2005, 105, 2999–3093.
79. S. Chen, R. Vijayaraghavan, D. R. MacFarlane and E. I. Izgorodina, J. Phys. Chem. B, 2013, 117, 3186–3197.
80. A. D. Buckingham, J. E. Del Bene and S. A. C. McDowell, Chem. Phys. Lett., 2008, 463, 1–10.

---