Unveiling the Borohydride Ion through Force-Field Development

The borohydride ion, BH4–, is an essential reducing agent in many technological processes, yet its full understanding has been elusive, because of at least two significant challenges. One challenge arises from its marginal stability in aqueous solutions outside of basic pH conditions, which considerably limits the experimental thermodynamic data. The other challenge comes from its unique and atypical hydration shell, stemming from the negative excess charge on its hydrogen atoms, which complicates the accurate modeling in classical atomistic simulations. In this study, we combine experimental and computer simulation techniques to devise a classical force field for NaBH4 and deepen our understanding of its characteristics. We report the first measurement of the ion’s activity coefficient and extrapolate it to neutral pH conditions. Given the difficulties in directly measuring its solvation free energies, owing to its instability, we resort to quantum chemistry calculations. This combined strategy allows us to derive a set of nonpolarizable force-field parameters for the borohydride ion for classical molecular dynamics simulations. The derived force field simultaneously captures the solvation free energy, the hydration structure, as well as the activity coefficient of NaBH4 salt across a broad concentration range. The obtained insights into the hydration shell of the BH4– ion are crucial for accurately modeling and understanding its interactions with other molecules, ions, materials, and interfaces.

Pitzer model on NaCl/NaOH mixture and extrapolation to pH=7:

Experimental data
We conducted osmolality measurements on ternary NaCl/NaOH solutions.This enabled us to verify the reliability of an extrapolation scheme that employs known properties of ternary NaCl/NaOH mixtures to predict the properties of pure binary NaCl solutions.The raw experimental data are presented in Tab.S1, in which the water activities and osmotic coefficients are also evaluated (Eq.1).

Advanced analysis
In order to develop a reliable empirical force field for BH − 4 , we have analyzed in detail the anion hydration by ab initio MD.We highlight the differences between the results from abinitio MD and a classical force-filed (FF) simulations, that would represent highly polar BH − 4 with ∆∆G solv ≈ −78 kJ/mol.
First, the radial distribution functions were determined for B-O w , B-H w , H-O w , and H-H w .Second, the directionality of the H-H w interaction (i.e., H-H 'hydrogen' bond) was determined via the angle-resolved g(r, α).Last, the 3D-density maps of water (i.e., of water oxygen and water hydrogen) around BH − 4 were calculated and compared.
From ab initio simulations, we have also analyzed the distribution of bond lengths and angles between BH − 4 and the hydration water molecules, and, importantly, also the partial charges in BH − 4 .We observed a charge transfer from BH − 4 to the neighboring water molecules, which further lowers the polarity of BH − 4 and the partial charges of its hydrogen atoms.
The radial distribution function g(r) was calculated and normalized as usual.The spatial distribution function g( r) is defined via Eq. 1, i.e., as a ratio of a local and a bulk density at a position r Bi .We note that the central particle is the boron atom (placed at [0,0,0]), with the best possible orientation of two hydrogen atoms, where the first one defines xaxis.With the second hydrogen, they define the reference xyz-coordinate system.In order to symmetrize and smoothen the 3D-spatial distribution function, we have averaged over all possible hydrogen atom permutations.We have calculated the spatial distribution of water oxygen and hydrogen with the spatial resolution of a bin = 0.25 Å, i.e., volume element In order to calculate the angle-resolved radial distribution function for water hydrogens, we proceeded in two steps (see Figure 5).First, the distance between boron and water hydrogen was calculated (r BHw ).Then the closest boron hydrogen to water hydrogen was found (r HHw ) and the angle H nearest -B-H w was calculated, see the sketch in Figure 5.The normalization is performed with respect to both, the angular (α) and spatial coordinate (r-distance).

Results
In this section, we compare the hydration structure of BH This difference could be roughly related to the signatures of typical force-field hydration of NH + 4 vs.CH 4 .In the former case, the solute-water H-bonding features are important.In the latter, the water H-bond network from bulk water is essentially preserved, just encapsulating the 'large' spherical solute (hydration of small hydrophobic solutes).
Comparing the 2D distribution functions in Figure S3, further clear differences between ab-initio and FF simulations are observed.First, the empirical force field exhibits much stronger orientational preferences (see yellow spots at low-α, r < 4 Å).Secondly, the B-H w (and consequently also H-H w ) distance is shorter by about 0.1 Å compared to the ab-initio result.See also the discussion related to the RDFs presented in Figure S2.
The strong orientation of water in MD FF and directionality of the H-H w bonds, which Table S3 presents the significantly different charge distribution within the BH − 4 anion pre-  dicted by ab-initio (MD) with charge transfer (q B = −0.56e,q H = −0.11e),which contrasts with (q B = +1.2e,q H = −0.55e)used in the FF simulation to match ∆∆G solv ≈ −78 kJ/mol (literature estimate).Table S3: Partial charges on the neutral borane (BH 3 ) and borohydride anion (BH − 4 ) as calculated at different levels of theory.The second column presents the results of gas-phase calculation, the third column presents the result from ab-initio MD of BH − 4 in a spherical water cluster.In the case of gas-phase, the partial charges are derived using the RESP method.In ab-initio MD simulation, the unrestricted (EPS) partial charges were calculated.Moreover, it was found that about 20 % of the charge (i.e., 0.2 e) was transferred to the water molecules in the first hydration layer, therefore also the rescaled values of the partial charges (total charge of BH − 4 q = −1e) are also shown (fourth column).The partial charges employed in the empirical FF of highly polar BH − 4 are shown in the last column for comparison.

QM: Hydration Free Energy
Hydration free energy from cluster continuum approach Table S4: Cluster-continuum approach (CCA) to calculate ∆G solv for different numbers (N wat = 0 -3) of explicitly involved water molecules (solvation energies in kJ•mol −1 ) calculated at MP2/aug-cc-pVTZ employing PCM with the Bondi radii.

Molecular dynamics (MD) simulations
The pair interaction potential between atoms V ij is modeled through the sum of Coulomb and Lennard-Jones (LJ) interactions where q i , q j are the charges of atoms i, j and r ij is the distance between these atoms.
We use the Lorentz-Berthelot combination rules to define the Lennard-Jones parameters where i, j correspond to the index of the atoms and ions.
We employ the SPC/E water model. 3The model assigns partial charges of −0.8476 and 0.4238 to oxygen and hydrogen, respectively.The water geometry is fixed at a bond length of 1.0 Å and a bond angle of 109.47 • using the LINCS algorithm. 4For Na + and Cl − , we use the Smith-Dang parameters. 5

Free energy of solvation
The solvation free energy of the BH − 4 ion is calculated using thermodynamic integration 6 where H λ is the Hamiltonian of the system, λ LJ and λ C are the LJ and charge transition coordinates, which are 0 in the initial state and 1 in the final state.The solvation path is split in two separate processes: first, a neutral van der Waals particle is created, which is assigned a charge in the second step.Along the transition path, the λ-dependent Hamiltonian is defined as We set the exponent k = 6 in this equation to avoid divergences.Integrations are performed through a 12-point Gaussian quadrature with λ ∈{0.00922, 0.04794, 0.11505, 0.20634, 0.31608, 0.43738, 0.56262, 0.68392, 0.79366, 0.88495, 0.95206, 0.99078}.For every value of λ, we perform a 400 ps simulation of which the first 50 ps are discarded for equilibration.
The ionic solvation free energy computed in the simulations is sensitive to the simulation scheme (system shape, periodic or finite system) and treatment of the electrostatic forces (Ewald sum, cut-off based, etc).Therefore for comparison with experimental data, several corrections have to be applied to the simulation data.The correction term accounting for finite system and ion size reads 7 where z is the ion valency and e the elementary charge.Here, R ion is the effective radius of the ion, which is estimated from the ion-oxygen radial distribution function, and ε r is the relative dielectric constant of the SPC/E water.The Wigner potential is ξ ew = −2.837279/L,where L is the simulation box size in nm.
Experimental values of the solvation free energy are usually given with respect to a hypothetical transfer of ions from the ideal gas phase of p 0 = 1 atm pressure to the ideal solution under pressure of p 1 = 24.6 atm, corresponding to a density of 1 mol/l.Thus, it is also necessary to include a correction term related to the compression of the gas ∆G press = N A k B T ln(p 1 /p 0 ) = 7.9 kJ/mol where N A and k B T are Avogadro's number and the thermal energy, respectively.
Hence, the total single-ion solvation free energy is given by ∆G solv = ∆G sim + ∆G fs + ∆G press . (5) Solvation free energy and diffusion coefficient The green symbols show the MD simulation results for solvation free energy.The partial charges and LJ interaction parameters of the borohydride ion are q B = 0.108 e, q HB = −0.277e, ε B = 0.4 kJ/mol, σ HB = 0.34 nm, and ε HB = 0.5 kJ/mol.As boron atoms are screened by hydrogen atoms, changing their radius does not affect the solvation free energy.

Self-diffusion coefficient of a single ion in water:
For the final LJ parameters of BH − 4 , the self-diffusion coefficient is calculated from an additional 100 ns NVT ensemble simulation of the single ions in different cubic box sizes (L = 3.0, 4.0, and 5.0 nm).All simulations are pre-equilibrated in the NPT ensemble before fixing the box size.The self-diffusion coefficient of the borohydride ion is calculated from the slope of the mean-square displacement by a linear fit as shown in Fig. S9.The diffusion coefficient corrected for system size effects is calculated using the formula: where L is the box length, D pbc the computed self-diffusion coefficient, D 0 the diffusion coefficient for infinite non-periodic systems, k B the Boltzmann constant, T the absolute temperature, η the solvent viscosity, and ξ ew = 2.837297 the self-term for the cubic lattice.
To correct for the low viscosity of the SPC/E water compared to the measured water viscosity, we report the scaled diffusion coefficients, 9

Figure S1 :
Figure S1: Non-ideality of NaBH 4 solution as quantified by (a) osmotic coefficient, (b) activity coefficient, and (c) activity coefficient derivative.Dashed lines represent the combined uncertainty arising from both experimental measurements and the fitting model.
− 4 anion obtained from ab-initio MD simulation and force-field that represents polar BH − 4 anion that reproduces the literature estimate of ∆∆G solv ≈ −78 kJ/mol.The radial distribution functions, presented in Figure S2, obtained from ab-initio and from force-field simulations are clearly different.Most importantly, at the ab-initio level the hydration water around BH − 4 is relatively dispersed, and the g(r) do not exhibit pronounced structure.Therefore, in ab-initio simulations, BH − 4 behaves similar to a large anion with a small charge density.In contrast, highly polar empirical force field exhibits sharp peaks in the first and second hydration shell, thus signaling very specific hydration of BH − 4 anion.

Figure S2 :
Figure S2: Radial distribution function between boron and water (left, B-H w black, B-O w red) and hydrogen atoms of BH − 4 and water (right, B-H w black, B-O w red).Hydration structures obtained from ab-initio MD simulations of a spherical domain (radius 7.5 Å described at QM level) embedded in TIP3P force-field water (top) and an empirical force-field (bottom, polar BH − 4 of ∆∆G solv ≈ −78 kJ/mol) MD simulation are compared.

Figure S3 :
Figure S3: Probability distribution of H w -H BH −4 orientations, defined by the α-angle, as explained in Figure8in the main text.The distance from a water hydrogen atom to the central boron atom is measured, and the angle is defined between the closest boron hydrogen, the boron atom, and the water hydrogen.The 2D-distribution obtained from ab-initio MD simulation shown is on the left, while the FF-based result (∆∆G solv ≈ −78 kJ/mol) is presented in the middle.Note that the visualization has been capped at a maximum probability of 10 for clarity, but the full peak height is displayed on the right.

Figure S4 :
Figure S4: Geometry of the first hydration layer of the BH − 4 anion obtained from ab-initio (left) and force-field (right, ∆∆G solv ≈ −78 kJ/mol) simulations.Note the significantly enhanced directionality of water-anion H-H w bonding in the case of the force field simulation, and bulk-like H-bonding structure in the case of the ab-initio simulation.

Figure S5 :
Figure S5: Spatial distribution functions of water hydrogen (white) and oxygen (red) atoms around the BH − 4 anion.The ab-initio results are shown on the top, and the force field simulation results (polar BH − 4 of ∆∆G solv ≈ −78 kJ/mol) at the bottom.The isocontour level of the spatial distribution is enhanced by 3× for O w and 2× for H w relative to the bulk density in the case of ab-initio MD, and 15× (both O w , H w ) in the case of FF simulations.

Figure S6 : 4 . 1
Figure S6: Microhydrated structures of Cl − , BH − 4 , and BF − 4 optimized on the MP2/augcc-pVTZ level of theory.While ordered and strongly directional interactions between water molecules and Cl − anion are observed, less structured hydration is observed for BH − 4 and BF − 4 .

Figure S7 :
Figure S7: Borohydride ion with a labeling notation of individual atoms.

Figure S8 :
Figure S8: The solvation free energy of BH −4 ion as a function of boron atom radius.The green symbols show the MD simulation results for solvation free energy.The partial charges and LJ interaction parameters of the borohydride ion are q B = 0.108 e, q HB = −0.277e, ε B = 0.4 kJ/mol, σ HB = 0.34 nm, and ε HB = 0.5 kJ/mol.As boron atoms are screened by hydrogen atoms, changing their radius does not affect the solvation free energy.

Figure S9 :
Figure S9: Mean square displacement BH − 4 ion using the SPC/E water model.Different color lines refer to different box sizes.Diffusion coefficients are obtained by a linear fit from 0.8 to 5 ns shown as dotted lines.

Table S1 :
Osmolality of ternary NaCl/NaOH mixed solutions as determined in VPO measurements at 310 K.The real ternary solutions, at finite NaOH concentrations, were employed during the fitting with the Pitzer model, while the binary NaCl solution was a target used to estimate the quality of the extrapolation.
We have measured the osmolality of ternary solution of NaBH4/NaOH at m NaOH = 1 mol/kg (pH ≈ 14).The raw experimental data are presented in Tab.S2, in which the water activity

Table S2 :
Osmolality of ternary NaBH 4 /NaOH mixed solutions as determined in VPO measurements at 310 K.The data served determining remaining (free) parameters in the Pitzer model.NaOH = 0.1 mol/kg were ignored since our experiments proved that the heat associated with slow NaBH 4 decomposition at pH ≈ 13 affected our VPO signal (hundreds of mOsm units).