Properties of Hydrogen Bonds in the Protic Ionic Liquid Ethylammonium Nitrate. DFT versus DFTB Molecular Dynamics

TL;DR

This study compares DFT and DFTB molecular dynamics simulations of ethylammonium nitrate, evaluating their accuracy in describing hydrogen bonds and infrared spectra, and finds DFTB suitable for such analyses.

Contribution

The paper assesses the performance of DFTB versus DFT in simulating hydrogen bonds in a protic ionic liquid, highlighting DFTB's strengths and limitations.

Findings

DFTB predicts similar average geometries to DFT.

Differences are notable in NH- and CH-stretching spectral regions.

DFTB fails to accurately describe nitrate vibrational fingerprints.

Abstract

Comparative molecular dynamics simulations of a hexamer cluster of the protic ionic liquid ethylammonium nitrate are performed using density functional theory (DFT) and density functional-based tight binding (DFTB) methods. The focus is on assessing the performance of the DFTB approach to describe the dynamics and infrared spectroscopic signatures of hydrogen bonding between the ions. Average geometries and geometric correlations are found to be rather similar. The same holds true for the far-infrared spectral region. Differences are more pronounced for the NH- and CH-stretching band, where DFTB predicts a broader intensity distribution. DFTB completely fails to describe the fingerprint range shaped by nitrate anion vibrations. Finally, charge fluctuations within the H-bonds are characterized yielding moderate dependencies on geometry. On the basis of these results, DFTB is recommend for the simulation of H-bond properties of this type of ionic liquids.

Figures

Tables

Table 1. Table 1: Comparison of position of the first peak (in Å) in the radial pair distribution function of the hexamer cluster as calculated with different methods. For comparison bulk neutron diffraction  Hayes11_3237 , DFTB  Addicoat14_4633 , force field MD  Song12_2801 , and Car-Parrinello (BP86) AIMD  Bodo13_144309 results are given as well. In the right-most column one finds the average HB angle.

	OH	ON	OC
BLYP	1.9	2.9	3.7	24.4
DFTB	1.9	2.8	3.7	29.7
Ref. Bodo13_144309	1.8	2.8	-	-
Ref. Hayes11_3237	2.4	3.0	-	-
Ref. Addicoat14_4633	-	2.5	3.0	-
Ref. Song12_2801	2.5	3.2	3.5	-

Table 2. Table 2: Average Mulliken charges along the DFTB and BLYP trajectories and partial charges of the force field (FF) from Ref.  Umebayashi08_125

	$N_{\mathrm{cation}}$	H	O
DFTB	-0.14e	0.28e	-0.73e
BLYP	-1.15e	0.53e	-0.01e
FF	-0.36e	0.31e	-0.635e

Equations

$$I(\omega )={\omega }^2⟨|\underset{0}{\overset{T}{\int }}\exp (-i\omega t)\vec{\mu }(t)\kappa (t)\mathrm{d}t{|}^2⟩.$$

$$q^A=Z_A-\int {\rho }_A(r)\mathrm{d}r.$$

$$N=\int \rho (r)\mathrm{d}r=\sum _{ab}{P}_{ab}{S}_{ab},$$

$$N_A=\int {\rho }_{\mathrm{A}}(r)\mathrm{d}r=\sum _{a\in A}{P}_{aa}+\frac{1}{2}\sum _{\begin{array}{c}a\in A\\ b\notin A\end{array}}{P}_{ab}{S}_{ab}.$$

Full text

††journal: Theoretical Chemistry Accounts

1 Introduction

Ionic liquids (ILs) are promissing candidates for novel applications, including the use as solvents and CO2 absorbers Firaha15_7805 ; Kerle09_12727 ; Welton99_2071 . Their special physico-chemical properties are determined by the nanoscopic structure and dynamics, which are governed mainly by strong Coulomb forces between the ions. For example, the ion pair lifetime can be linked directly to transport properties as shown by Zhang and Maginn Zhang15_700 . There also exist directional hydrogen bond (HB) interactions between ions of opposite charge within a wide range of HB strengths hunt15_1257 . Although being much weaker compared to Coulomb interactions, they impact macroscopic properties decisively. A prominent example is the reduction of the viscosity by the presence of HBs as shown in Ref. Fumino08_8731 .

Experimentally, HB interactions are accessible by far infrared (FIR) spectroscopy, as was discussed in detail by Fumino et al. Fumino12_6236 , as an extra blue-shifted band if compared with the effect of dispersion interaction. Furthermore, in the mid-infrared (MIR) region the H-bonded NH- or CH-stretching vibrations are red-shifted with respect to their non-H-bonded counterparts Roth12_105026 . The signatures of H-bonding and Fermi resonance interactions can be investigated in detail using nonlinear spectroscopic techniques, which are capable to unravel different phase and energy relaxation rates nibbering07_619 ; Zhuang09_3750 ; Chatzipapadopoulos15_2519 .

Numerical simulations can be used to support the various spectroscopic evidences for the existence of HB interactions. Thereby, one can distinguish between cluster and bulk liquid setups. Although these systems are rather different at first glance, gas phase spectroscopy of ILs in dependence on cluster size has revealed local motifs that give rise to bulk like spectra Johnson13_224305 . Structural similarities have also been found in simulations and were assigned to the dominance of strong short-range Coulomb interactions Bodo13_144309 . Cluster models are particularly well-suited for a normal mode analysis, corresponding to a minimum structure at zero Kelvin. The resulting frequencies are commonly scaled by a method-dependent factor and the spectra are artificially broadened to fit the experimental line shapes. Still, the pronounced anharmonicity of HB dynamics might require to calculate explicit potential energy surfaces giese06_211 . In addition, there exist trajectory-based normal modes approaches that can account for temperature effects mathias11_2028 , which also have been applied to ILs Thomas14_024510 . However, the standard approach especially for bulk simulations makes use of linear response theory, i.e. the IR spectrum is calculated by Fourier transformation of the dipole autocorrelation function may11 . The latter is obtained by sampling MD trajectories from a canonical ensemble ivanov13_10270 . Of course, the quality of the spectra depends strongly on the underlying potential energy surface driving the dynamics.

Classical MD simulations use many body interactions parametrized via force fields. There are many applications to structural and dynamical properties of ILs Maginn09_373101 ; DelPopolo04_1744 , discovering, for example, mesoscopic segregation behavior Bernardes14_6885 ; Shimizu14_567 . Concerning IR spectroscopy and in particular signatures of HB dynamics the use of force fields is problematic, since the parametrization of the latter is often targeted to thermodynamic quantities, see, e.g., Ref. zentel16_234504 . Here, DFT-based ab initio molecular dynamics (AIMD) provides the proper frame for more accurate simulations Marx10_ ; Bodo13_144309 ; kirchner15_202 . For instance, the structure and HB dynamics of methylammonium nitrate was investigated in detail by Zahn et al. using Car-Parinello MD Zahn10_124506 . Employing radial, angle and spatial distribution functions, preferred ion orientations and HB properties were investigated and the behavior of ion caging was observed. The latter has been discussed as a general feature of IL dynamics Schroder11_024502 . For small system sizes it is in principle possible to obtain AIMD trajectories long enough to calculate the power spectra Wendler12_1570 through velocity autocorrelation functions or, in combination with a charge localization scheme, IR spectra via dipole autocorrelation functions Thomas15_3207 .

Density functional-based tight-binding (DFTB) is a computationally less demanding alternative, thus allowing for larger system sizes Elstner98_7260 ; Yang07_10861 . DFTB doesn’t require empirical input and provides self-consistent Mulliken charges, which account for polarization effects. It has been successfully applied to a large class of problems, including biological systems Elstner06_316 , as well as to study structural properties of ILs zentel16_234504 ; Addicoat14_4633 .

In this contribution we address the performance of DFTB with respect to its description of HB dynamics and IR spectra. This is done for the test system ethylammonium nitrate (EAN), [(C2H5)NH${}_{3}]^{+}$[NO${}_{3}]^{-}$, which is a prototypical protic ionic liquid (PIL). The reliability of the results will be judged against data obtained by DFT/BLYP AIMD simulations. To facilitate the AIMD calculations the system is restricted to a hexamer cluster of six ion pairs. Previously, DFT and force field calculations have been compared for EAN clusters of various sizes by Bodo et al. Bodo12_13878 ; Bodo13_144309 . These authors put a special emphasis on HB geometries of energy optimized structures. In DFT structures they found asymmetric HB networks with only two out of three possible HB formed. In contrast force field simulations yielded more symmetric structures and saturated HBs. Furthermore, they compared a harmonic normal mode spectrum against Raman measurements and concluded on signals from NH-stretching vibrations involved in HBs of different strength. The force field itself had been parametrized by Song et al. by combining various parameter sets with additional ab initio data Song12_2801 . Radial distribution functions from MD simulations were found in reasonable agreement with high-energy X-ray diffraction data. Interestingly, OH-distances and HB angles in force field MD simulations are substantially larger than the DFT ones from energy optimized structures Bodo12_13878 and from later AIM simulations Bodo13_144309 .

DFTB has been compared to force field simulations for the related PIL triethylammonium nitrate in our previous work zentel16_234504 . Preceding bulk simulations, geometric HB correlations for a hexamer cluster have been compared for DFTB and DFT trajectories. This gave support for the reliability of the DFTB method, which then was used to quantify geometric correlations and to obtain IR spectra for the liquid phase. In particular for the FIR range, where experimental data are available, DFTB yielded excellent results whereas the force field failed to reproduce the experiments.

Alkylammonium cations provide the means to tune the network of HBs formed between cations and anions via the number of alkyl chains fumino09_3184 ; zentel16_ . Whereas triethylammonium nitrate is capable of forming a single HB only, the present system, EAN, in principle facilitates the formation of three HBs. In the following the HB geometries and IR signatures of EAN will be investigated using cluster MD simulations. Thereby, DFTB and DFT trajectory data will be compared, which will lend support for DFTB as an efficient yet accurate method for simulation of the HB network in PILs such as EAN. The paper is organized as follows: In the next section a brief overview on the used theoretical methods and computational setups is given. Results are presented in Section 3, focussing on geometries, IR spectra, and charge distributions. Finally, a conclusion is provided in Section 4.

2 Theoretical Methods

2.1 Molecular Dynamics

A cluster consisting of 90 atoms (six ion pairs) with similar initial structure as the hexamer cluster studied by Bodo et al. Bodo12_13878 , is used for the following DFTB and AIMD simulations. In a first step the cluster has been geometry-optimized using DFT with the BLYP exchange correlation functional and the 6-31+G* basis set. Further, Grimme’s D3 dispersion correction has been employed Grimme10_154104 . Starting from the minimum structure, Langevin dynamics was performed with the target temperature being 300 K and a Langevin damping time of 1 ps. The stochastic perturbations introduced by the Langevin thermostat are small enough to yield spectra equivalent to the ones from the standard approach (microcanonical trajectories sampled from an canonical ensemble) Gottwald12_ . The total NVT trajectory length has been 60 ps and the time step 0.4 fs. For calculating the IR spectrum, it is split into trajectories of 10 ps length. The DFT AIMD simulations, in the following called BLYP, were performed with Terachem, version 1.5K Ufimtsev09_2616 .

Starting from the same initial BLYP geometry, DFTB simulations were done using third-order density expansion Gaus11_931 and the 3-ob Slater-Koster parameter set Gaus13_338 . The DFTB simulations are carried out using the dftb+ software, Version 1.2 Aradi07_5678 . First, a canonical ensemble was simulated for 5 ps employng the Nosé-Hoover thermostat with target temperature 300 K. This was followed by an energy optimization. Subsequently, the optimized structure is used to initiate a canonical ensemble run for 100 ps. From this trajectory starting geometries for seven microcanonical (NVE) trajectories were sampled randomly. The microcanonical trajectories had a length of 10 ps and a time step of 0.4 fs.

DFT and DFTB trajectories were analyzed with respect to the presence of HBs. Here, a HB was defined between the nitrogen-bound hydrogen and its closest oxygen atom, without a distance cutoff. This way each of the six cations were assigned three HB partner oxygens from the anions. The meaningfulness of this choice is analyzed below.

2.2 IR Spectra

The IR spectrum, $I(\omega)$, is obtained by Fourier transformation of the dipole moment trajectory, →μ(t), using a Kaiser window function, $\kappa(t)$, with Kaiser parameter 10, i.e.

[TABLE]

Here, the average is taken with respect to the canonical ensemble. The window function ensures that the dipole fluctuations are damped to zero at the integration boundaries. The Kaiser parameter controls the time scale of the decay and is chosen large enough to ensure that the spectra are not artificially broadened. All spectra plotted are the moving averages of 21 steps to smooth the curves. In BLYP the total cluster dipole is calculated directly from the wave function, whereas in DFTB the Mulliken charges are used to calculate the total dipole trajectories of the cluster. Note that the dipole in an ionic systems is not uniquely defined, as discussed in Ref. Brehm12_5030 . In the context of IR spectroscopy the center of mass of the cluster is a plausible reference, as this point is fixed for all vibrational normal modes Thomas15_3207 . In addition to the total cluster dipole, the contribution of the $\mathrm{NH_{3}}$ group atoms to the total dipole are extracted from DFTB Mulliken charges.

For comparison vibrational normal modes are calculated for both methods. In case of DFT the Turbomole software package has been used Ahlrichs89_165 ; turbomole6.5 ; Deglmann04_103 .

2.3 Population Analysis

As compared to standard force fields DFTB has the advantage of delivering atom-centered point charges (Mulliken charges), which are not only determined self-consistently, but they are polarizable in the sense that they adopt to the actual environment of the atom zentel16_234504 .

In general the partial charge of an atom $A$ with atomic number $Z_{A}$ is defined by

[TABLE]

Depending on the method used the partial charges can vary because the atomic densities, $\rho_{A}(r)$, can be defined in different ways. A simple approach to obtain partial charges is the Mulliken population analysis. Here, the total electron density, $\rho(r)$, is written in terms of the density matrix elements, $P_{ab}$, and the atomic basis functions, which after integration over space gives the total number of electrons

[TABLE]

where $S_{ab}$ is overlap matrix element and the indices $a$ and $b$ refer to atomic basis functions. The Mulliken population $N_{A}$ of an atom $A$ is defined as

[TABLE]

To gain insight into the charge distribution within the HB, Mulliken analyses are performed for the respective N, O, and H atoms using a set of methods, including DFTB, DFT, and coupled cluster theory (CCSD(T)); Turbomole is used for the CCSD(T) calculations Ahlrichs89_165 ; turbomole6.5 ; Haser89_104 ; Hatting00_5154 .

3 Results

3.1 Cluster Structure and HB Geometries

The radial pair distribution functions, $g_{\mathrm{ON}}(r)$, $g_{\mathrm{OH}}(r)$, and $g_{\mathrm{OC}}(r)$ have been calculated from the distances between the anionic O to the cationic N, the H bound to it, and the alkyl C atoms, respectively. DFTB and BLYP results are compared in Fig. 2. The convergence with respect to the trajectory length was checked for the DFTB simulations. Note that due to the cluster structure, $g(r)$ approaches zero for larger distances. In general these pair correlation functions are rather similar for DFTB and DFT clusters. The largest difference between DFTB and DFT is visible for $g_{\mathrm{OH}}(r)$ around 2.5 $\mathrm{\AA}$, which is due to slightly different orientations of the $\mathrm{NO}_{3}^{-}$ molecules with respect to the N-H bond, bringing a second O atom closer to the N-H bond (for values of first peak, see Tab. 1).

HBs can be detected as distinct peaks at small distances between two atoms involved in HB. Indeed, $g_{\mathrm{OH}}$ and $g_{\mathrm{ON}}$ show maxima at short distances, but there is no such contribution in $g_{\mathrm{OC}}$. Thus HBs only develop via the cationic nitrogen atom to the anion, and not via the alkyl chains.

Cluster results are available from the work of Bodo * et al.*. In Ref. Bodo12_13878 they reported distances obtained for an equilibrium (geometry optimized) hexamer cluster structure. Here, average OH bond length of 1.86 Å and ON distances of 2.9 Å(DFT/$\omega$B97X-D) have been found. Car-Parrinello AIMD simulations of the same hexamer Bodo13_144309 yielded rather similar values, i.e. 1.8 Å and 2.8 Å for the average OH and ON distance, respectively. These values actually agree very well with those of bulk phase simulations reported in the same reference. This has been attributed to the similarity of the first solvation shell due to the prevalence of strong Coulomb interactions. There are other bulk phase investigations as well. For instance, an analysis of neutron diffraction data has been reported in Ref. Hayes11_3237 ; cf. Tab. 1. The ON distance compares rather well with the present findings. The OH distance is smaller by 0.5 Å, both in DFTB and BLYP. Comparing with force field MD simulation data of Song et al. Song12_2801 this could be a bulk effect, although such a conclusion would be at variance with the AIMD results of Ref. Bodo13_144309 . In this respect we note that the bulk DFTB results reported in Ref. Addicoat14_4633 are markedly different from experimental and MD data. To conclude, the present DFTB setup yields average HB geometries, which are consistent with both BLYP and Car-Parrinello results.

As can be seen in Fig. 1 the cations orient themselves such that the ethyl chains are on the outer border of the cluster and the polar head group is close to the anions and forms HB to oxygens from different anions. The HBs were identified from the trajectory based on the shortest $r_{\mathrm{OH}}$ found, thus all nitrogen bound hydrogens are supplied with a H-bonding partner oxygen atom. To verify this procedure, all HB pairs are analyzed as well by the heavy atom distance, $r_{\mathrm{NO}}$, combined with an angle criterion. Applying the classifications of moderate HB strength of $r_{\mathrm{NO}}<3.2~{}\mathrm{\AA}$ and $\alpha<50^{\circ}$ from Ref. hunt15_1257 , 15.6 % and 10.1 % of the configurations are not HBs in DFTB and BLYP, respectively. Interestingly, in both cases in less than 1 % of the samples distance and angle criteria are violated at the same time. Moreover, in 29% (DFTB) or 25 % (BLYP) of the cases one has twofold HBs, where one oxygen is acceptor of two H atoms. These twofold HBs do not show increased OH distances or HB angles, compared to the other HBs, suggesting that they are of similar strength. Thus we conclude that in the present cluster the NH groups are mostly saturated and form an asymmetric network where an O atom can partake in two HBs, leaving some O atoms without H-bonding partners.

These findings are in contrast to the DFT minimum energy structure reported in Ref. Bodo12_13878 . Although in that case also an asymmetric HB network has been observed, where only two of the possible three HBs were formed, but there were no twofold HBs. This difference could be due to the finite temperature used in the present simulation.

In the following the N-H$\ldots$O HBs are discussed, putting emphasis on geometrical correlations. To that end, HB coordinates $r_{1}=0.5(r_{\mathrm{NH}}-r_{\mathrm{OH}})$ and $r_{2}=r_{\mathrm{NH}}+r_{\mathrm{OH}}$ are introduced and plotted against each other in Fig. 3 (left panel) for both methods. If the distances $r_{\mathrm{NH}}$ and $r_{\mathrm{OH}}$ were fully uncorrelated, the results would lie on a straight line with the slope -2, plotted for reference in grey color. Within a HB the H atom is pulled towards the acceptor atom. Thus, $r_{\mathrm{OH}}$ is reduced as $r_{\mathrm{NH}}$ increases. This results in deviations from the linear relation between $r_{1}$ and $r_{2}$ and shifts the curve at constant $r_{1}$ to larger $r_{2}$ values. Furthermore, the closer $r_{1}$ is to zero, the stronger is the bond, as the hydrogen is located exactly between the two heavy atoms, corresponds to the strongest HB possible limbach06_193 .

To compare observed correlations to other H-bonding systems, the valence bond model of Pauling is used. Here, starting from the HB distances, $r_{\mathrm{NH}}$ and $r_{\mathrm{OH}}$, bond orders are defined as $p_{i}=\exp(-(r_{i}-r_{i}^{\mathrm{eq}})/b_{i})$ with $i=\mathrm{\{NH,OH\}}$ pauling47_542 . Under the constraint that the sum of the two bond orders $p_{i}$ must be equal to one, the two coordinates depend on each other and the proton/H atom transfer can be described by a single coordinate Limbach04_5195 . This path is drawn in Fig. 3 as a blue line. The equilibrium distances, $r_{i}^{\mathrm{eq}}$, are obtained from single molecule gas phase geometry optimizations. The value for the N-H bond is taken as the average of the three ammonia N-H bond lengths of the optimized structure, which are $r_{\mathrm{NH}}^{\mathrm{eq},\mathrm{DFTB}}=1.029~{}\mathrm{\AA}$ and $r_{\mathrm{NH}}^{\mathrm{eq},\mathrm{BLYP}}=1.036~{}\mathrm{\AA}$. The O-H bond length is $r_{\mathrm{OH}}^{\mathrm{BLYP}}=0.991~{}\mathrm{\AA}$ in BLYP and $r_{\mathrm{OH}}^{\mathrm{DFTB}}=0.982~{}\mathrm{\AA}$ in DFTB. The bond order decay parameters were changed simultaneously starting from the values reported for HBs in crystal structures steiner98_7041 until a good visual fit with the DFTB data was observed, resulting in $b_{\mathrm{OH}}=0.3~{}\mathrm{\AA}$ and $b_{\mathrm{NH}}=0.33~{}\mathrm{\AA}$. These values were not optimized separately for BLYP, as the DFTB results describe the BLYP correlations equally well. Overall, the HB distance correlations shown in Fig. 3 (left panel) are strikingly similar in the two methods. In passing we note that a similar analysis was performed for trietylammonium nitrate using DFTB bulk MD simulations, where the resulting bond order decay parameters were found to be $b_{\mathrm{NH}}=0.355~{}\mathrm{\AA}$ and $b_{\mathrm{OH}}=0.321~{}\mathrm{\AA}$ zentel16_234504 . Interestingly, judging from the geometric correlations, the HBs in the EAN cluster are weaker than in the above bulk structure.

Finally, we comment on the linearity of the HBs in the EAN cluster. Fig. 3 (right panel) shows the distribution of O atom positions in the plane defined by the N, H, and O atom making the HB. Although being rather similar the deviations from linearity are more pronounced for DFTB as compared with BLYP. However, on average DFTB reproduces the EAN structures from BLYP simulations very well as can be seen from the summary provided in Tab. 1.

3.2 IR Spectra

The IR spectra calculated from BLYP and DFTB trajectories using Eq. (1) are shown in Fig. 4. In panel (a) the complete spectral range is given, which reveals a qualitative agreement. That is, most peaks are present in both spectra, but the intensities vary strongly. This discrepancy is particularly striking in the 750 $\mathrm{cm^{-1}}$ to 1900 $\mathrm{cm^{-1}}$ region. To what extend the spectra deviate from one another due to differences in intensity only, or due to a different underlying nuclear dynamics, has been scrutinized using a normal mode analysis. In Fig. 5 DFTB and BLYP mode frequencies are correlated based on the overlap between the respective normal mode displacement vectors. Leaving aside the overlap, one finds a reasonable linear correlation for the range up to 900 cm*-1* and for the N-H and C-H stretching range, i.e. here the density of states is rather similar for the two methods. Strong deviations are observed for the fingerprint region. Taking into account the overlap we notice the following: In the region up to about 350 cm*-1* the overlap is in general below 50 %, from 400 to about 900 cm*-1* the overlap is considerably larger, in some case reaching almost 100 %. In the region of the stretching vibrations (see inset of Fig. 5), there is a reasonable overlap, but only the three NH-stretching modes with highest frequency reach about 100 %. The problematic fingerprint region shows a good overlap for those modes, which match in frequency. In contrast, those modes which differ in frequency also have a small overlap. In the following we will compare the spectrum in the different spectral regions separately.

The FIR region shown in panel (b) is shaped by intramolecular vibrations as well as intermolecular dispersion and H-bonding interactions. The experimental spectrum for bulk EAN from Ref. fumino09_3184 is also given. It exhibits a double peak structure with a broad maximum at around 200 $\mathrm{cm^{-1}}$ and a smaller peak around 420 $\mathrm{cm^{-1}}$. Using DFT/B3LYP cluster calculations, Fumino et al. fumino09_3184 assigned the low-frequency band to asymmetric and symmetric HB stretching vibrations at 197 and 128 $\mathrm{cm^{-1}}$, respectively, and the high-frequency band to a N-C-C deformation vibration of the cation. The shoulder around 250 $\mathrm{cm^{-1}}$ was found to be due to ethyl-nitrogen torsional motions. Note that this assignment was independent on the cluster size (from one to six pairs) and has been confirmed by the present BLYP normal mode calculations. In view of the correlation in Fig. 5 it is interesting to note that the spectrum as such is rather similar even though the overlap between the normal mode vectors is only modest.

Comparing theory and experiment, both DFTB and BLYP reproduce the observed double peak structure. The main peak is blue-shifted by about 25 $\mathrm{cm^{-1}}$ and its shoulder is not clearly resolved. The small differences between theory and experiment could be due to the method or due to different H-bonding in the cluster as compared with the bulk liquid. The high-frequency peaks are broader than in the experiment, with DFTB predicting again a blue shift of about 25 $\mathrm{cm^{-1}}$. As far as the different widths are concerned one should note that in the gas phase cluster, the ethyl chains are oriented to the outer part of the structure and are therefore able to explore the phase space more freely than in a bulk system, which could lead to a broadening of the peak.

H-bonding also manifests itself in red-shifted NH-stretching vibrational frequencies of those modes involved in a HB with respect the free NH case. Details of the corresponding spectral region are shown in Fig. 4c. The IR spectrum obtained from the DFTB trajectory shows a broad plateau-like feature from 2900 to 3100 $\mathrm{cm^{-1}}$, while the DFT signal shows a double peak structure with a maximum at 3042 $\mathrm{cm^{-1}}$ and a distinct side band at 3166 $\mathrm{cm^{-1}}$. From the Mulliken charges of the DFTB simulation it is possible to extract just the dipole from a specific subgroup of atoms to better assign spectral features. The resulting signal from the NH3 group dipole is included in Fig. 4. The signal has two peaks at 2837 $\mathrm{cm^{-1}}$, and 2940 $\mathrm{cm^{-1}}$, and is very broad, an indication of HBs with various strengths. The remaining intensity in this range should be due to CH-stretching vibrations.

In order to judge the predictions in this spectral range, we compare with the Raman spectra of liquid EAN published in Ref. Bodo12_13878 and analyzed based on hexamer DFT/$\omega$B97X-D normal mode calculations. In liquid EAN NH-stretching vibrations are found across the whole absorption band, i.e. from 2900-3300 $\mathrm{cm^{-1}}$, reflecting the distribution of HB strengths. In addition, the cluster normal modes contain signatures of very strong HBs (around 2800 $\mathrm{cm^{-1}}$) and free NH-stretching vibrations (around 3400 $\mathrm{cm^{-1}}$). CH- stretching vibrations are assigned to the range between 2900-3000 $\mathrm{cm^{-1}}$). The observation that the low-energy part of the spectrum is dominated by H-bonded NH-stretching vibrations agrees with the present DFTB findings. Clearly, the finite size cluster cannot capture the broad distribution of HB strengths present in the liquid. In the BLYP case the absorption band is narrower hinting at less structural flexibility of the cluster (see also angular distribution in Fig. 3). Further, the first moment of the BLYP band appears to be blue-shifted as compared to DFTB. Given the fact that the highest frequency NH-stretching modes in the correlation plot, Fig. 5, coincide, this evidences a more pronounced anharmonicity in DFTB as compared with BLYP.

Most interesting, however, is the range between 750 $\mathrm{cm^{-1}}$ and 1900 $\mathrm{cm^{-1}}$, which shows rather pronounced differences in Fig. 4(a). According to the Raman analysis in Ref. Bodo12_13878 in this range one expects the NO${}_{3}^{-}$ bending ($\sim$720 $\mathrm{cm^{-1}}$), symmetric stretching ($\sim$1045 $\mathrm{cm^{-1}}$), and asymmetric stretching ($\sim$1400 $\mathrm{cm^{-1}}$) vibrations. First, we assign this region of the IR spectrum according to the BLYP normal mode analysis as follows (in parentheses the position and intensity of the strongest peak) : 685 - 720 cm*-1* – NO${}_{3}^{-}$ bending (717 cm*-1*, 14 km/mole), 780 - 795 cm*-1* – NO${}_{3}^{-}$ out-of-plane bending (780 cm*-1*, 32 km/mole), 998 - 1032 cm*-1* – NO${}_{3}^{-}$ symmetric stretching (998 cm*-1*, 32 km/mole), 1489 cm*-1* – NO${}_{3}^{-}$ asymmetric stretching (1489 cm*-1*, 413 km/mole), 1561 - 1680 cm*-1* – NH3 bending (1654 cm*-1*, 152 km/mole). The range from 1002 - 1472 cm*-1* contains rather mixed type of vibrations, involving mostly CH3 and CH2 bending motions and NO${}_{3}^{-}$ asymmetric stretching vibrations. The largest intensity is at 1432 cm*-1*(463 km/mole). Taking into account that the AIMD results contain the effects of anharmonicity and finite temperature, one could tentatively assign the spectrum in Fig. 4(a) as follows: 785 cm*-1* – NO${}_{3}^{-}$ out-of-plane bending, 1000 cm*-1* – NO${}_{3}^{-}$ symmetric stretching, double peak around 1320 cm*-1*– CH3 and CH2 bendings and NO${}_{3}^{-}$ asymmetric stretching vibrations, double peak around 1580 cm*-1* – NH3 bending vibrations. Compared to the Raman data there are some deviations, which are difficult to discuss since the experiment has been performed in the liquid phase. Decomposing the power spectrum of a Car-Parrinello simulation on the same size cluster Bodo et al. Bodo13_144309 obtained the NO${}_{3}^{-}$ out-of-plane bending at 730 cm*-1*, the NO${}_{3}^{-}$ symmetric stretching at 962 cm*-1*, and the NO${}_{3}^{-}$ asymmetric stretching yielded a broad feature in the range 1280 - 1370 cm*-1*. Overall, one can conclude that the present BLYP simulations are consistent with available experimental and simulation data.

Apparently, the DFTB spectrum looks different in this spectral region and the question arises whether this is an effect of intensity, vibrational frequency or both. In order to scrutinize this issue, we have modified the correlation plot between DFTB and BLYP frequencies, taking into account only points where the BLYP intensity of the associated transition is above a certain threshold. If the threshold is such that all peaks in the fingerprint region are accounted for, the figure doesn’t look much different from Fig.  5(a). The result shown in panel (b) of Fig.  5 has been obtained for a threshold of 100 km/mole, which sets the focus on the strongest peak in the fingerprint region. From these observations we conclude that for the lower intensity peaks the vibrational frequencies match reasonably between BLYP and DFTB, but the intensities are different. For the region of the most intense band, which is due to CH3 and CH2 bending motions and NO${}_{3}^{-}$ asymmetric stretching vibrations, both peak positions and intensities are different. In view of the favorable comparison of BLYP with previous simulation and experimental results, it must be concluded that DFTB does not yield a correct IR spectrum for the nitrate anion’s vibrations. In passing we note that this finding is in accord with the conclusions drawn in Ref. otte07_5751 , which reported deviation between DFT and DFTB up to 10% for NO-related vibrations of nitric acid.

3.3 Charge Distributions

The IR spectral regions of HB-related frequencies are well reproduced using DFTB. This suggestes to use the DFTB-Mulliken charges to characterize HB-related atomic charge fluctuations in the thermal ensemble. First, we inspect the Mulliken charges of those atoms involved in H-bonding given in Tab. 2. They have been averaged over the whole trajectory and all H-bonded pairs. The DFTB values appear to be closer to the force field parameters from Ref. Umebayashi08_125 than to the BLYP ones. Closer inspection reveals that this is due to the combination of a small basis set and using Mulliken charges. Because a sufficiently large basis set is computationally not affordable for the BLYP simulations, we will consider only DFTB trajectory data in the following.

In particular we will be interested in the distributions of Mulliken charges on the N, O, and H atom of all HBs in the cluster with respect to the internal HB coordinates $r_{1}$ and $r_{2}$. In order to estimate the accuracy of DFTB in this respect, gas phase calculations of Mulliken charges for a single ion pair have been performed using different methods, i.e. DFTB, BLYP, and CCSD(T) (for basis sets, see figure). The resulting charges upon HB elongation are shown in Fig. 6. Taking CCSD(T) as a reference we notice that DFTB reasonably reproduces the CCSD(T) behaviour in case of the N atom, i.e. sign and slope match. For the O atom the sign matches as does the slope for short HB distances. The strongest deviation in slope is seen for the H atom. As far as BLYP is concerned, it is interesting to note that it fails to give the correct sign of the charge at the N atom.

In Fig. 7 we present the changes of Mulliken charges with respect to the average values sampled from the DFTB trajectory. The following trends can be observed: The charge at the O atom shows a broad distribution with respect to both coordinates, but no obvious correlation. In contrast, for the H and N atom the distributions are narrower and at the same time a coordinate dependence is clearly visible. The H atom’s charge increases (i.e. less electron density) towards smaller absolut $r_{1}$ and smaller $r_{2}$, values, and at the same time the charge of the N atom decreases. As far as the H atom is concerned it can be anticipated, based on Fig. 6 that the actual change should be even more pronounced.

4 Conclusion

Molecular Dynamics simulations of a hexamer cluster of the protic ionic liquid EAN were performed with the goal being a critical assessment of the validity of the DFTB approach. To this end a comparison with DFT/BLYP as well as CCSD(T) results has been made. Averaged HB geometries agreed quantitatively for DFTB and BLYP and are in line with earlier experimental and theoretical publications. With respect to the linearity of the HB BLYP predicts a slightly less flexible structure. HB coordinate correlations are rather similar for DFTB and BLYP, both being in the range of weak correlations. The IR spectra calculated from the cluster dipole trajectories show considerable differences between the two methods. In the far IR region the spectral features due to H-bonding are similar and close to available bulk experimental data. In the stretching region, DFTB yields a somewhat broader spectrum pointing to a larger flexibility of the cluster. Most critical is the fingerprint region, where DFTB clearly fails to give reasonable frequencies and intensities for the NO${}_{3}^{-}$-related vibrations. A comparison of the normal modes between the two methods confirms this conclusion.

To sum up, DFTB can be considered as computationally less demanding alternative as compared to DFT if the focus is set on the properties and dynamics of HBs.

Acknowledgements.

The authors thank the Deutsche Forschungsgemeinschaft (DFG) for financial support through the SFB 652 and A. Wulf for supplying the experimental FIR spectra.