Molecular dynamics simulations of ammonium/phosphonium-based protic ionic liquids: influence of alkyl to aryl group †

The variation of the center atom in the cation from an N to a P-atom leads to improved physiochemical properties of protic ionic liquids (PILs) which are suitable for electrolyte applications. We present an atomistic simulations study to compare the eﬀect of an alkyl or aryl group on trioctylammonium triflate ([HN(Oct) 3 ][TFO]) and triphenylammonium triflate ([HN(Ph) 3 ][TFO]) with their phosphonium analogues. We have computed the binding energy from quantum chemical calculations and physical properties such as the viscosity and the electrical conductivity of PILs from molecular dynamics simulations. The influence of the aromatic character in PILs is found to be significant to the physical properties. Gas phase quantum chemical calculations on clusters of ion pairs have revealed the presence of C–H/ p interactions in aromatic PILs along with hydrogen bonding. The variation in strength of the ion-pair affinities is examined using electric-current correlation and velocity autocorrelation functions. The qualitative differences observed are due to the aromatic rings and change in the central atom of the quaternary cation from an N to a P-atom, substantiated quantitatively by diffusion coefficients and electrical conductivities. The relatively weaker ion-pair interactions and low binding energy ( (cid:2) 73.34 kcal mol (cid:2) 1 ) lead to the highest electrical conductivity in [HP(Ph) 3 ][TFO].


Introduction
Electrochemical properties at wide temperature ranges popularized ionic liquids (ILs) as successful electrolytes in batteries, dye-sensitized solar cells, capacitors and fuel cells. [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16] In the past few decades, tremendous effort has been put towards the investigation of the compatibility of ILs with different electrochemical devices in anhydrous or aqueous phases. 4,9,17 Parallel to the experimental findings, substantial support through a molecular-level understanding has been provided over the years 18,19 to explore ionic mobility, vibrational signatures of various interactions [20][21][22] and ion hopping. 23,24 Recently, the existence of one or more than one solid-solid phase transition in ionic liquid plastic crystals has drawn attention due to the specific contribution of various solid phases to the electrolytic properties. [24][25][26][27][28][29][30][31] Protic ionic liquids (PILs) were favoured for their potential use in portable and long life cycle energy devices over acid electrolytes mainly due to a high energy density, a large electrochemical window and a wide range of thermal stabilities. In the search for proton conducting ILs, the applicability of ammoniumbased PILs has been well examined in the past decade. Rana et al. 32 synthesized tri-butylphosphonium cation-based ammonium analogue PILs and demonstrated their facile ionic conductivity compared to ammonium-based PILs. The use of phosphonium-based ILs as an electrolyte is interesting due to their competent/improved physical and electrochemical properties with the ammonium-based PILs. 32, 33 Luo et al. 34 explored a comparative study of tri-octylphosphonium triflate and triphenylphosphonium triflate with ammonium analogues. The authors found a higher cation conductivity in phosphonium salts compared to their ammonium analogues which was attributed to weaker hydrogen bonding in phosphonium cation-based ILs. Apart from these protic ILs, the effect of the variation of the center atom species of the cation from an N to a P atom on the physicochemical properties of aprotic quaternary cation-based ILs also revealed similar results. For example, tri-n-ethylpentylphosphonium bis-trifluoromethanesulfonyl amide ([TEPP][TFSA]) showed a lower viscosity and higher self-diffusion coefficient compared to its ammonium analogues (i.e. tri-n-ethylpentylammonium bis-trifluoromethanesulfonyl amide ([TEPA][TFSA])). 35 Recently, Salem et al. 36 Fig. 1). The computational details of the gas phase calculations and MD simulations are provided in the next section.

A. Gas phase
The initial structures of the ion pairs of these PILs were constructed using GaussView 41 software. We have carried out geometry optimizations on five different initial configurations (for the monomer, dimer, trimer, and tetramer respectively) which were chosen randomly from the MD simulations' trajectories. For e.g., five different initial configurations of the [HN(Ph) 3 ][TFO] dimer are shown in Fig. S1 of the ESI. † The geometry optimized structures correspond to the minimum energy among the five different configurations analysed.
(i) Calculations using atom-centered basis sets. For the geometry optimization of an isolated ion pair (monomer) in the gas phase, a M06/aug-cc-pVDZ level of theory was implemented using the Gaussian09 code. 42 Vibrational frequencies were calculated at the same level of theory and the absence of any imaginary frequency confirmed that the ion pairs were in their true minimum energy configuration.
(ii) Calculations using plane-wave basis sets. DFT calculations for clusters of two (dimer), three (trimer) and four ion pairs (tetramer) for all of these PILs were carried out using CP2K software. 43 The effect of the core electrons and nuclei on the valence electrons were considered by using Geodecker-Teter-Hutter (GTH) pseudopotentials. 44 Triple-z double-polarized basis sets with an energy cutoff of 280 Ry were employed to treat the valence electrons. A Perdew, Burke, and Ernzerhof (PBE) exchangecorrelation functional was employed. 45 The van der Waals interactions were incorporated through Grimme's DFT-D3 scheme. 46 A cluster of these PILs was quenched in an isolated condition with convergence criteria of 10 À5 and 10 À7 for nuclear coordinates and gradient of wave functions, respectively. The minimized structures were used for further analysis.

B. Liquid phase
Classical molecular dynamics simulations were performed on the liquid phases of these PILs using the LAMMPS package. 47 We used 256 ion-pairs for each system and the simulations were carried out at 393 K, 417 K, 441 K and 465 K. All bonded interaction parameters were taken from the OPLS-AA force field. 48 Non-bonded potential parameters for these salts were taken from the work of different groups. In our previous study, for alkylammonium-based PILs we used the van der Waals parameters of Chang et al. 49  , OPLS-AA Lennard-Jones parameters were used. In order to include the effect of electronic polarization, which plays a crucial role in determining the structure and dynamics of ILs, 51-53 the atomic site charges for these PILs were calculated using the DDEC/c3 program. 54,55 For this purpose, an isolated ion pair of these salts was optimized within the DFT framework using the CP2K package. 43 The minimized coordinates were used to generate the valence electron density and were saved in a cube file. This density information was used as an input to the DDEC/c3 code to obtain atomic charges. Other details of the calculations are the same as those described in the DFT calculations for clusters. The atomic site charges obtained for these salts are tabulated in Tables S1-S3 of the ESI. † The remaining nonbonded and bonded force field parameters are provided in Tables S4 and S5 of the ESI. † The total ion charges were in a range between AE0.70 to AE0.78e which can be understood as an effect of the electronic polarization and charge transfer present in these systems. Fig. 1 The chemical structures of tri-octyl/tri-phenyl ammonium triflate and tri-octyl/tri-phenyl phosphonium triflate PILs.
Energy minimization was performed on the initial configuration created using Packmol. 56 A distance cutoff of 12 Å was employed to calculate the pairwise interactions in real space. Long-range interactions were computed through a particleparticle particle mesh (PPPM) solver with a precision of 10 À5 . All of the C-H covalent bonds were constrained using the SHAKE algorithm as implemented in LAMMPS. 47 The velocity Verlet algorithm was used to integrate the equations of motion with a timestep of 1 fs. The Nosé-Hoover thermostat and barostat 57,58 were employed to control the temperature and pressure of the system. All the systems were equilibrated in an NPT ensemble for 15 ns. The computed densities from the NPT simulations are given in Table 1. A trajectory of 48 ns was generated for each system in the NVT ensemble to calculate several equilibrium properties. The systems were visualized in VMD. 59 3 Results and discussion

Binding energy and structure
The binding energy (BE) is estimated as follows where E IL , E Cation , and E Anion represent the energy of an isolated ion-pair, cation and anion, respectively. The binding energies are corrected for the basis set superposition error (BSSE These hydrogen bonding interactions are further characterized using the radial distribution functions (RDFs) obtained from MD simulations of 256 ion-pairs of ILs. RDFs between the acidic hydrogen atom attached to the cationic N-atom or P-atom with anionic oxygen atoms are shown in Fig. 3 , respectively. It is due to the influence of charge transmission from the N/P to the electron withdrawing phenyl ring compared to positive electrostatic induction by the octyl chain. These results clearly support a previous finding by Fumino et al. 20 where an enhancement in the electron releasing nature of the ammonium cation of protic molten salts/PILs leads to weaker hydrogen bonds. A similar trend in aryl to alkyl substitution in these PILs is seen on N-H stretching frequencies from the experimental findings of Luo et al. 34 due to a weakening of the hydrogen bond. Moreover, these RDF peak profiles are in excellent agreement with hydrogen-bond distances obtained from quantum chemical calculations. The high peak intensity at a relatively lower distance of N-HÁ Á ÁO interaction in [HN(Ph) 3 ][TFO] among all PILs leads to a strong ion-pair binding affinity. Similar observations are seen in the center of mass RDFs calculated for cation-cation and cationanion interactions (Fig. S2 of ESI †).
In order to examine the distribution of anions around the cation, the spatial density distribution maps of anions (S-atom and O-atoms) around the center of mass of the cation were calculated from MD simulations at 393 K. Fig. 4 shows the spatial density maps where the anion distribution over the acidic hydrogen is more dense for aliphatic ILs compared to aromatic PILs. Moreover, the presence of bulky octyl chains leads to a more condensed anion distribution above or below the acidic hydrogen (N-H or P-H bond vector). A random orientation of the The interactions of anions with the ring hydrogen atoms of the cation are insignificant. The influence of the above ion-pair interactions and binding energy on physico-chemical properties is illustrated in the next section by calculating the shear viscosity and dynamics of PILs.

Shear viscosity
The shear viscosity of these systems is computed using the equilibrium Green-Kubo relation, 61 where V is the system volume, and P is the symmetric, traceless part of the pressure tensor. Pressure tensors were stored at each step for a 4 ns duration and twelve such independent trajectories were generated. The block average of these trajectories was used to compute the shear viscosity from the time integral  of the stress autocorrelation function (eqn (2)). The viscosity (in mPa s) calculated from the MD simulations of tri-octyl/triphenyl ammonium triflate and tri-octyl/tri-phenyl phosphonium triflate PILs is shown in Table 2

Electric current autocorrelation function.
The electrical conductivity s GK is calculated using the electric current autocorrelation function which can be defined 62,63 as where, the j(t) is the electric-current function equal to P N i¼1 q i v i ðtÞ; q i and v i represent the charge and velocity of atom i at time t].
The electric current autocorrelation function is the sum of the cation-cation autocorrelation ( j ++ ), anion-anion autocorrelation ( j ÀÀ ) and cation-anion cross correlation function ( j +À ). 64 These individual contributions are computed and shown in Fig. 5. The well depth seen in the decay of the electric-current function  (see inset of Fig. 5) for all PILs indicates the caging of ion-pairs. The relative difference in the ion-pair cage strength is further examined by computing the center of mass velocity autocorrelation function.

Velocity autocorrelation function.
In order to characterize the dynamical processes in PILs, the normalized center of mass velocity autocorrelation function (VACF) for cations and anions is calculated (see Fig. 6). The first minimum is found to be deeper for the tri-aryl cations compared to tri-alkyl cations (see inset of Fig. 6a for [HN(Ph) 3

][TFO] and [HP(Ph) 3 ][TFO]).
A significant decrease in well depth with an increase in the bulky nature of quaternary ammonium/phosphonium cation is also seen in our recent work. 60 The observation from cation VACF in the well depth variations suggests that the rattling motions seem more significant in [HN Fig. 6a. The exponential decay of the anion VACF is found to be different in tri-aryl cations than tri-alkyl cations (Fig. 6b).
In order to quantitatively verify these observations, the diffusion coefficients of the ions (D A in Â10 À6 cm 2 s À1 ) are calculated. The D A calculated from MD simulations for tri-octyl/ tri-phenyl ammonium triflate and for their phosphonium analogues is provided in Table S7 of the ESI. † For all PILs, the D A of anions in the respective PIL is higher than compared to cations. For all temperatures, [TFO] anions diffuse faster than the cations which could be due to the bulkier nature of the cations compared with the anions. Moreover, the ionic conductivity within the approximation of independent ion motion can be calculated from diffusion coefficients using the Nernst-Einstein relation. 63 The s NE for all the ILs is given in Table S8 of the ESI. † We have calculated the electrical conductivity (s GK ) using the Green-Kubo relation at various temperatures (see Table 3 among all the PIL which shows the least intermolecular ion-pair interactions and a low ion-pair binding affinity. Moreover, an increase in temperature leads to an increase in s GK . From 393 K  , the s GK increases by a factor of 10, 11 and 14 respectively. Thus, an increase in temperature from 393 K to 465 K gives a tenfold or more than tenfold enhancement in s GK . Moreover, the ratio of s GK and s NE is calculated (see values in parenthesis of Table 3) which illustrates the ion-pair association. Except 393 K, the ratio of s GK and s NE suggests that ion-pair association in ammonium-based PILs is higher compared to their phosphonium analogues. It is due to more significant hydrogen-bonding interactions seen in ammonium-based PILs. These comparative trends seen in ionpair association of ammonium salts and their phosphonium analogues are in excellent agreement with experiments 34 on these PILs.

Conclusions
The quantum chemical calculations and MD simulations provide a molecular understanding of various properties of tri-alkyl/ tri-aryl ammonium/phosphonium-based PILs. The difference in distances obtained from quantum chemical calculations for hydrogen bonding and C-H/p interactions clearly distinguishes the strength of ion-pair interactions for respective ammoniumbased PILs versus their phosphonium analogues. The RDFs for N-HÁ Á ÁO and P-HÁ Á ÁO interactions reveal that hydrogen bonding interactions are weaker for phosphonium-based PILs. Similar observations are seen from the spatial distribution of anions (O/S atoms) over the acidic hydrogen of the cation. The distribution of anions around the acidic hydrogen of the cation is found to be the least for [HP(Ph) 3  The dynamics of the ion-pairs are characterized by calculating the electric current autocorrelation function. The ion-pair association behavior is seen from the ratio of s GK and s NE which suggests that ion-pair association is higher in ammonium-based PILs compared to their phosphonium analogues. These results provide a molecular understanding of the superior electrochemical nature of phosphonium-based PILs over their ammonium analogues. This study can be further extended to investigate proton diffusion through ab initio MD methods and to explore the role of P-atom versus N-atom in proton diffusion.

Conflicts of interest
There are no conflicts to declare.