УДК 541.135 FORCE FIELD OF TETRAFLUOROBORATE ANION FOR MOLECULAR DYNAMICS SIMULATION: A NEW APPROACH

Design of new electrical energy storage devices including supercapacitors as well as an optimization of existing ones require not only new electrolytes, but also the deep and complete understanding of the processes occurring in the electrolyte solutions. Spectral techniques and classical molecular dynamics simulation (MDS) have gained a reputation as a reliable tool for such tasks. The starting point of any MDS is a choice or development of the force fields for all simulated particles. The combination of vibrational spectroscopy and molecular dynamics technique can provide a thorough understanding of the structure and dynamics of the ionic subsystem. In this connection, the reproduction of the vibrational spectra should be added to the requirements for the force fields of the most common electrolyte components. Many modern supercapacitors are based on organic electrolytes consisting of non-aqueous aprotic solvents such as acetonitrile, propylene carbonate and γ-butyrolactone and quaternary ammonium salts with tetrafluoroborate and hexafluorophosphate as anions. The purpose of the current work is to develop a new force field for tetrafluoroborate anion (BF4) able to reproduce not only translational diffusion in acetonitrile medium, but also the spectral properties of this ion in a condensed phase. Since found in the literature force fields of BF4, cannot satisfy these requirements, there were performed intensive quantum chemical calculations of BF4 at the M06-2X/6-311++G(d,p) level of theory to construct the potential energy surface with respect to the B-F bonds and F-B-F angles followed by evaluating corresponding intramolecular potential constants. Combining the obtained bond and angle force constants with partial charges on B and F atoms calculated at the same level of theory, and literature values of Lennard-Jones parameters, a new force field model for BF4 anion was created. Based on the carried out MD simulations of the BF4 ion in an infinitely diluted acetonitrile solution, it was proved that the obtained resulting model is capable to reproduce both transport and intra-ion vibrational properties of the tetrafluoroborate anion.


Introduction
Rapid technology increase, necessity of finding an alternative to non-renewable resources stimulated the development of new devices for storing electric energy. Supercapacitors are one of such kind devices [1][2]. Supercapacitors have several advantages such as big power density, fast charging / discharging time, long lifetime etc. In the same time they have several ways to improve its technical characteristics first of all capacitance and service voltage. Capacity and power density depends on containing electrolyte [1], therefore deep understanding of undergoing processes in electrolyte medium and interaction mechanism between solution components are the keys to successful development of new, more effective devices.
Spectral methods have proven to be a reliable tool for studying ion-molecular systems [3][4][5][6], but, like many experimental methods, for a comprehensive understanding obtained result, they need to involve computational methods, such as classical molecular dynamics simulation (MDS) [7][8]. Therefore, the development of electrolyte component force fields capable of accurately reproducing spectral properties is an important task of modern electrochemistry. Being a component of a number of electrolytes for supercapacitors [1], tetrafluoroborate anion (BF 4 -) can be called one of the most priority candidates for developing such a force field. A significant part of the manufactured supercapacitors offer devices with organic solvents, mainly with propylene carbonate or acetonitrile [1]. High demand for acetonitrile (AN) as a supercapacitor component is a reason of choice this substance as a solvent of investigated systems. https://doi.org/10.26565/2220-637X-2019-  Based on the foregoing, the goal of this work was to create a model of the tetrafluoroborate anion force field, able to reproduces the spectral and dynamic properties of this ion in acetonitrile.
To achieve this goal, first we reviewed existing approaches to modeling tetrafluoroborate anion as a component of an electrolyte solution (Section 1). Since classical MD simulation requires an information about the equilibrium structure and atomic partial charges for any multiatomic particle, the corresponding parameters were obtained by quantum chemical calculations as well as the vibrational spectrum of BF 4 in vacuum (section 2). In the Section 3 the details of the MD simulations are described along with the results of simulation of BF 4 vibrational spectra for the most common models in the literature. Section 4 presents the results of the potential energy surface (PES) scanning of B-F bonds and the angles of F-B-F followed by evaluation of the corresponding intramolecular potential constants. Finally, in the section 5 the validation of the developed force field of tetrafluoroborate anion is presented based on the anion diffusion coefficient in the infinitely diluted acetonitrile solution.
CG anion of tetrafluoroborate is one particle with a total mass of five atoms and the corresponding charge. The advantage of such models is a significant saving computational resources, which allows to compute larger systems and longer phase trajectories. On the other hand, within the framework of CG approach a detailed study of ionic subsystem structure on atomic level is impossible.
Polarizable models allow one to take into account various effects of electrostatic interaction, but the existence of several ways [12] of implementation does not give complete confidence in the use of this approach.
Rigid full-atomic models save computational resources due to fixed bond lengths and angle values and provide "atomic resolution" of the ionic subsystem study. In the same time, it is impossible to obtain vibrational spectra of tetrafluoroborate anion with restricted intraionic structure.
Based on the foregoing, within current investigation, we decided to use full-atom soft force field models of tetrafluoroborate as an approach able to completely solve the problems posed to us. As a result of our literature review, models ABS [25], Liu / Wu [23][24] and CL&P [26] identified as the most popular tetrafluoroborate force fields.

Quantum-chemical calculations of the optimal geometry, partial atomic charges, and the vibrational spectrum of BF 4 in a vacuum
All quantum chemical calculations were carried out using the Gaussian 09 [27] software package were performed at the M062X/6-311++G(d,p) level. The specific charges on the atoms were calculated by the Breneman ChelpG method [28].
As a result of quantum-chemical calculations of the tetrafluoroborate anion in vacuum, following equilibrium geometry parameters were obtained: the B-F bond length equal to 1.4092033966 Å, and the F-B-F angle equal to 109.47122063°. The correct tetrahedral form of the ion was set as the initial requirement for the results of quantum-chemical calculations. For boron and fluorine atoms the partial charges of +1.188568 and -0.5471420.4569 were evaluated.
The obtained vibrational modes, their expected intensities in the IR spectrum, as well as a description of the oscillations are given in table 1.
The data presented in table 1 are in good agreement with the literature [6].

Vibrational spectra for the most common models in the literature
Based on conducted literature review, next three force field models of the tetrafluoroborate anion were selected for spectral properties efficiency testing: 1) model of Andrade (Jones de Andrade), Boes (Elvis S. Böes) and Stassen (Hubert Stassen) [25], hereinafter referred to as ABS; 2) the Liu / Wu model, first proposed by Liu [23], and then refined by Wu in [24]; 3) the model of Lopez (José N. Canongia Lopes) and Padua (Agílio A. H. Pádua) [26] hereinafter referred to as CL&P. The ABS model is one of the very first soft full-atom tetrafluoroborate models. The purpose of discussed research was to create force fields for modeling of room temperature ionic liquids (RTIL) based on alkyl-imidazoline cations with tetrachloroaluminate and tetrafluoroborate anions as counterions. The AMBER [29] force field was used as the basis of the model, but since the AMBER force field is developed for modeling biological objects [29], it has no parameters of atoms untypical for biological field. The Lennard-Jones parameters of boron atom were taken from the more universal DREIDING force field [30], but parameters for fluorine atoms were taken AMBER [29]. Quantum chemical calculations were performed by GAMESS [31] program on the UHF / 6-31G (d) level of the theory; boron and fluorine atoms partial charges were set +0.8275 and -0.4569, respectively. Work of Andrade [25] was the first attempt to obtain interionic constants for B-F bonds and angles F-B-F of tetrafluoroborate anion using PES. This procedure was performed using the AMBER utilities. The ABS force field was used in a number of works [18][19]22], and was even used as the basis for another model of the force field of the tetrafluoroborate anion -Liu / Wu.
The Liu / Wu model was also used in a number of studies [17,20,24] was first presented in [23], to improve the force fields of imidazolium based ionic liquids. The force field in the discussed article is similar to the tetrafluoroborate model of Andrade [25], namely the same forms of intramolecular and intermolecular potentials as well as AMBER [29] force field as source were used. Comparing to ABS model [25], the Liu / Wu model has different BF bond equilibrium value (it was taken from crystallographic data [32]), different level of quantum-chemical calculations (HF/6-31+G(d)), and therefore different atomic partial charge values (B +1.1504; F -0.5376). Lennard-Jones parameters, B-F bond potential constants, equilibrium value and potential constant of F-B-F were taken from ABS model [25] with rounding. Despite the author's attention mostly was attracted to the cations, the advantage of Liu's work [23] is presence of calculated vibrational spectra (using classical MD-simulation) and their comparison to experimental ones by vibrational frequencies.
Wu's work [24] is devoted to improving the developed Liu's force field [23] and investigation of various concentrations ionic liquids acetonitrile solutions. Compared to [23], the anion's force field parameters are represented with greater accuracy and closer to [25] values.
The third discussed force field of tetrafluoroborate, CL&P, is well-known for ionic liquid force field development. This force field is a combination of the OPLS-AA [33] and AMBER [29] fields with further author's modifications.
Despite the fact that in a number of works [34][35][36] the authors themselves directly declare the tetrafluoroborate anion as a component of the simulated system and provide links to publications devoted to the development of the ionic liquids force fields [26,[37][38][39], the data for tetrafluoroborate force field parameters were found in supplementary information for [26] as DL_POLY [40] force field file. The values of σ and ε are equal to σ and ε from ABS [25] with rounding to the fourth decimal place, the partial charges on the boron and fluorine atoms equal +0.96 and -0.49, respectively. As in all previously considered models, the F-B-F angle equilibrium value is 109.5° and B-F bond equilibrium value is 1.394 Å, but CL&P intermolecular potential constant values of bonds and angles are 1.3 times and 1.5 times larger respectively.
In current work we have tested the force field models described above with respect to their ability to reproduce spectral properties of discussed anion in classical molecular dynamics simulation.
Intermolecular interactions were calculated as a sum of short-range (U SR ) and long-range (U Сoulomb ) interactions: The particle-mesh Ewald (PME) summation was used to calculate the long-range electrostatic interactions [41].
The Lennard-Jones potential (2) were used to calculate short-range interactions. It can be written both in classical and A-B form: 12 6 where ε ij is the depth of the energy well, and σ ij is the distance of zero interaction energy between the two particles, ij Intramolecular potentials for bonds and angles were calculated according to formulas (3) and (4), respectively: where K r is the intramolecular potential constant, r 0 is the bond length equilibrium value.
where K θ is the intramolecular potential constant, θ 0 is the angle equilibrium value. 1-4 interactions were disabled. Intermolecular and intramolecular potential parameters of discussed force field models (ABS, Liu / Wu and CL&P) are given in tables 2 and 3, respectively.
The vibrational spectra were calculated in vacuum and in a solvent medium by Fourier transform (5) of the autocorrelation functions of linear atomic velocities where vv C is normalized autocorrelation function of the linear velocity of the particle center-of-mass. All classical molecular dynamics simulations in a solvent medium were performed in NPT ensemble. Velocity rescaling thermostat [44] was used to maintain a constant temperature with T = 298.15K and the time constant value of 0.1 ps. The pressure was kept equal to of 1 atm using a Berendsen barostat [45] with the time constant of 10 ps. The simulated system consisted of one BF 4 ion and 500 AN molecules. Acetonitrile force field model of Koverga [46] was chosen for all MD simulations. Each simulated systems underwent through the following stages: energy minimization, thermodynamic equilibration (1 ns) and calculating properties (1 ns).
All molecular dynamics simulations in vacuum were performed in NVT ensemble with the volume of the simulated cubic cell equal to 100 nm 3 .
Obtained vibrational spectra for the ABS, Liu / Wu, and CL&P models are shown in Figures 1, 2, and 3, respectively.  Based on the presented results (see Fig. 1-3), for all models, vibrational modes positions in vacuum and in solvent are almost the same. Placing the ion into solvent medium simplifies the vibrational spectrum, namely three distinct pics in vacuum convert into one broad peak in acetonitrile accompa-nied by the shift of its position of ~ 5-10 cm -1 towards the long wave length range comparing to the central one in vacuum. This behavior occurs for all solvent spectra of all calculated models. Therefore for saving computational resources all further calculations of vibrational spectra by MD simulation were carried out in vacuum taking into account that the positions of triple peaks in the range of 600-700 cm -1 were represented by one middle peak. On presented spectra (see Fig. 1-3) there are oscillations in the low-frequency region (0-100cm -1 ) corresponding to inhibited translations and librations. The presence of such oscillations proves correct methodology of the molecular dynamics simulation.
The good convergence of the ABS and Liu / Wu models vibration frequencies is explained by the similarity of the intramolecular potential constant values: Liu / Wu values [23][24] were taken from [25]. This fact confirms the primary influence of intramolecular potentials over intermolecular in influence on vibrational spectra.
Both figures 1 and 2 have four modes: around 310 cm -1 , 360 cm -1 , 600 cm -1 and 1200 cm -1 ., respectively. These vibrational modes can be assigned to the modes from the table 1: 347.19 cm -1 , 519.06 cm -1 , 778.86 cm -1 and 1105.67 cm -1 , respectively. Based on given values, all modes of the Liu/Wu and ABS models show lower frequencies comparing to results of quantum chemical calculations for all modes except the last one (1200 cm -1 ). Frequency values obtained for CL&P force field model ( figure  3) are closer to table 1 values: around 400 cm -1 , 450 cm -1 , 700 cm -1 and 1200 cm -1 ., respectively. Comparing to Liu/Wu and ABS models the CL&P molecular dynamics spectrum has the same value of the last frequency mode (the forth one) but higher and therefore closer to quantum spectrum (see table 1) values.
Summarizing the results mentioned above, one can conclude that the ability of literature force field models [23][24][25][26] to reproduce vibrational spectrum of tetrafluoroborate by classical molecular dynamics simulation is quite questionable. Despite CL&P model is more effective than Liu/Wu and ABS, none of chosen force field are able to fully reproduce ion spectrum.

Scanning the surface of the potential binding energy of energy of the B-F bond and F-B-F angle
As it was shown previously, intramolecular potentials make the main contribution to the vibrational spectra, and none of discussed literature force field models can fully reproduce the vibrational spectra of tetrafluoroborate, it was decided to scan the potential energy of B-F bonds and F-B-F angles to obtain intrinsic molecular potentials of their own. Three methods of scanning potential energy were selected: 1) scanning with additional optimization of geometry (hereinafter, designate them as PES-1); 2) scanning without additional optimization of geometry (hereinafter, designate them as PES-2); 3) scanning without further optimization over four B-F bonds and two vertical F-B-F angles (hereinafter, we denote them as PES-3). All the quantum chemical calculations were carried out using the same level of theory as described in the Section 2.
The last method (PES-3) was used due to two reasons. First, it was an attempt to solve the problems with a drop-down point during the bond / angle scanning for highly symmetric molecules. Second, based on the table 1 the vibration III is a symmetric vibration of all four B-F bonds and there are no single or asymmetric vibrations of the B-F bonds. In this connection, for a better representation of vibration III, it is necessary to analyze the simultaneous symmetric vibration of all four B-F bonds. Figures 4 and 5 show the results of applying PES-1, PES-2, PES-3 to scan the energy profile for the B-F bond and the F-B-F angle, respectively, as well as the results of approximating the obtained dependences by the formulas (3) and (4). (3) and (4) describes well obtained PES-1, PES-2, PES-3 dependences. All three methods show similar function behavior in the region near to the optimum and geometric parameters values close to the corresponding optimal geometry as well.

Figures 4 and 5 demonstrates that equations
The values obtained during PES-1, PES-2 and PES-3 were approximated by the next equation: where c is a constant guaranteeing a zero value of the potential at the optimum point; K is K r or K θ ; x is r or θ; x 0 is r 0 or θ 0 for (3) or (4), respectively. The constants values obtained during the approximation, as well as their errors, are given in tables 5 and 6.  As expected the values of C constants are negligible (tables 5 and 6), that indicates the methodological correctness of the applied procedure.
Obtained sets of the intramolecular potential parameters (tables 5 and 6) were used to test their effectiveness in reproducing the vibrational spectrum of the tertafluoroborate anion in MD simulation. Liu/Wu force field model was chosen for this purpose. The calculation was carried out in vacuum according to the procedure described above. The results of the MD simulation of the vibrational spectra of BF 4 anion with the intramolecular potentials based on PES-1, PES-2 and PES-3 approaches are presented in Figure 6.  6. Molecular dynamics vibrational spectra of the BF 4 anion with intramolecular potentials based on PES-1 (red), PES-2 (blue) and PES-3 (green) and the frequencies of vibrational modes (grey vertical lines) from Table 1.
Observation of the figure 6 allows one to conclude that the PES-3 shows the best result for oscillation III (Table 1). PES-2 and PES-3 describe well oscillation II, and the results of PES-1 are equidistant from oscillations I and II. The task of reproducing vibration IV was not posed, since the tetrahedral form of BF 4 is preserved for this vibration. Because none of the used methods can provide the intramolecular potentials parameters able to completely reproducing of the vibrational spectrum of the BF 4 anion, it was decided to use obtained values (see tables 5 and 6) as a starting point for adjusting procedure allowing to evaluate the intramolecular parameters able to show the best agreement between the vibrational spectra from molecular dynamics simulations and quantum chemical calculations. The adjusting procedure consisted of a multidimensional analysis of the intramolecular potential to determine the values of K r and K θ , which give the best possible pick locations. The values of partial atomic charges and Lennard-Jones potentials were not used in this analysis because of the insignificant effect of intermolecular interactions on the vibrational spectrum of BF 4 as compare with intramolecular potentials (see section 3). Despite the fact that the equilibrium values of bonds and angles belong to intramolecular parameters, they were also not used as tunable parameters because of the importance of keeping the equilibrium geometry of the anion equal to that obtained in quantum calculations (see section 1). As a result, the values of K r = 410000 kJꞏmol -1 nm -2 , K θ = 710 kJꞏmol -1 ꞏrad -2 were obtained. All other parameters of the force field for of BF 4 anion were adopted as following. Equilibrium geometry parameters of the tetrafluoroborate anion (r(B-F) = 0.14092033966 nm and θ(F-B-F) = 109.47122063°) were calculated in section 2, as well as partial atomic charges (q(B) = +1.88568, q(F) = -0.547142). Lennard Jones parameters were taken from Liu/Wu model [23][24] (σ(B) = 0.35814 nm, σ(F) = 0.31181 nm; ε(B) = 0.3975 kJꞏmol -1 , ε(F) = 0.2552 kJꞏmol -1 .
The vibrational spectrum of BF 4 in vacuum obtained by MD simulation based on the final force field is shown in Figure 7.  Figure 7 demonstrates that the developed force field allows to reproduce the vibrational mode III as well as to achieve the equidistant position of the vibration frequencies for modes I and II.

Validation of the developed force field by the diffusion coefficient of the tetrafluoroborate ion
Force field model of ion as a component of electrolyte solution requires not only good representation of vibrational spectra, but also good representation of dynamic properties, such as the diffusion coefficient.
The procedure of diffusion coefficient calculation for BF 4 was identical to the solvent's one. The simulated system for this simulations included 1 ion and 500 solvent molecules. Figure   The results of calculating the diffusion coefficient of BF 4 according to the Einstein equation and the Green-Kubo formula (7) are shown in Table 7. As it follows form the table 7, the average calculated value of BF 4 diffusion coefficient in infinitely diluted acetonitrile is 2.30ꞏ10 -9 m 2 ꞏs -1 . Taking into account the correction factor of 1.17 for solvent diffusion, the obtained diffusion coefficient of the BF 4 anion can be estimated as 2.69ꞏ10 -9 m 2 ꞏs -1 . An experimental value of an ion diffusion coefficient can calculated from the limiting ionic conductivity of an ion by using the Nernst-Einstein equation where z i is the charge of an ion, Λ 0 m , is the molar limiting ionic conductivity of an on at infinite dilution.
Taking the limiting ionic conductivity of the BF 4 anion equal to 107.48 Smꞏcm 2 ꞏmol -1 [56], the value of BF 4 diffusion coefficient was estimated as 2.87ꞏ10 -9 m 2 ꞏs -1 . This means that developed force field of BF 4 is able to reproduce experimental diffusion coefficient of the anion in AN medium within the relative error of 6.3%.

Conclusions
Development of new supercapacitors and improving the existing ones requires thorough investigation of processes occurring in an electrolyte solution. Combination of spectral methods with molecular dynamics simulation can become a reliable tool for solving this kind of tasks.
This work is devoted to the development of BF 4 force field able to reproduce not only the dynamic, but also the spectral properties of the ion in acetonitrile media.
As a thorough literature review of BF 4 force field models, three models were chosen for further research: Liu/Wu [23][24], ABS [25] and CL&P [26]. Discussed force fields for classical molecular dynamics simulations of the anion in both vacuum and solvent medium. There is no significant difference between BF 4 spectra in vacuum and solvent medium was revealed, which indicates the prevalence of intra-molecular potentials over inter-molecular potentials in their influence on the vibrational spectrum. Since all discussed force fields are not able to reproduce spectral propertied of the ion with required accuracy, the potential energy surface scanning (PES) of bonds B-F and angles F-B-F were applied by using quantum chemical calculations at the M062X/6-311++G(d,p) level of theory.
Three different techniques of PES generation gave three different sets of intra-molecular parameters. To check these parameters, additional quantum chemical calculations of BF 4 at the M062X/6-311++G(d,p) level of theory were carried out to evaluate the atomic partial charges and equilibrium geometry. Then, taking into account the Lennard Jones potential parameters from Liu/Wu [23][24] model we have tested new sets of the intramolecular parameters in classical molecular dynamics simulation. Despite all obtained spectra of PES-generated force fields are not able to fully reproduce spectrum of BF 4 -, revealed tendencies of the intramolecular potentials influences on specific bands of the vibrational spectrum. The intermolecular parameters obtained from the PES analysis were used as starting point for their next generation. In addition, revealed relationship between the changes of the intramolecular constant values and the shift of the corresponding vibrational frequencies were used as a guide in multidimensional analysis of the intramolecular parameters. This allowed us to develop a new force field for the BF 4 anion with the best possible spectrum reproduction within the framework of classical molecular dynamics modeling.
The final BF 4 force field model were developed as a force field model of the electrolyte solution component, therefore reproduction of ionic diffusion properties is as important as reproduction of vibrational spectrum. For this purpose, the difference between diffusion coefficient of Koverga acetonitrile [46] and experimental value was calculated, and further this information were took during the validation of BF 4 diffusion. The relative error of the final BF 4 force field diffusion coefficient is 6.3%.