A transferable force-field for alkali metal nitrates

1.1.1. Lithium nitrate

1.1.2. Sodium nitrate

1.1.3. Potassium nitrate

1.1.4. Rubidium and caesium nitrate

Crystallise in a complex trigonal structure (Z = 9; P31, a = 10.55 Å, c = 7.47 Å (RbNO₃ at 298 K) [23] and a = 10.902 Å, c = 7.74 Å (CsNO₃ at 296 K)) [24].

Standard thermodynamic data is quoted in the literature for some of the polymorphs. Rao and Rao [25] discuss the enthalpies of transformation that are available. They conclude that the most reliable data are when there is agreement between the differential thermal analysis and calorimetric values. There are two examples known: for KNO₃ (phase II (Pnma) → phase I ($R\bar{3}m$)); ${\rm{\Delta }}H=5.0$ kJ mol⁻¹) and for RbNO₃ (phase IV (P31) → phase III ($Pm\bar{3}m$)); ${\rm{\Delta }}H=4.0$ kJ mol⁻¹). These small differences make reproducing the order of stability of the polymorphs a serious challenge. Figure 1 (below) shows the structures for the four space-groups where all the sites are fully occupied that have been used in the fitting process.

Zoom In Zoom Out Reset image size

The GULP [29, 30] code was used to fit the force-fields using the relaxed fitting algorithm [31]. Details of the code, including the methods used to calculate the long-range electrostatic summations, can be found in [30]. An iterative procedure was used to obtain the charges on the nitrogen and oxygen atoms. First the individual nitrates were fitted to experimental data to obtain interactions between the alkali metal ions and the atoms of the nitrate anion. Then the nitrogen charge, q_N , and short-range oxygen-oxygen intermolecular interactions were refitted to data for the observed low temperature polymorphs of all the alkali earth nitrates together, assuming that the intramolecular interactions of the nitrate ions and the intermolecular interactions (including electrostatic interactions) between the nitrate ions were the same for all the alkali metal nitrates. The oxygen charge, q_O , is fixed by the requirement that the total charge of the nitrate anion be −1. Then the alkali metal - nitrate interactions were re-fitted for each compound. This process was iterated until a self-consistent solution was obtained.

A considerable amount of experimental data is available for fitting the force-field. This includes structural data for many of the polymorphs, often at several temperatures (see the earlier discussion); elastic properties including the thermoelastic constants (d[log C_ij ]/dT) [32]; dielectric data [33] (more details are given in the supplementary information). Although partial occupancy can be modelled in GULP using a fractional occupation number, structures with partial occupancy were not used in the fitting process. Corrections need to be made to the available data before they can be compared with simulations. Lattice parameters obtained from experiment will be at a finite temperature and contain zero-point motion effects. Both these must be corrected for to compare with a lattice statics simulation (most often used in fitting). Lattice energies can be obtained from experiment and compared with the simulations, but individual enthalpies of formation are not easily usable since they require a model of the elements in their standard states.

We have performed density functional calculations using the CASTEP [34] (Cambridge Serial Total Energy Package) program to calculate the structures and properties of selected phases of the nitrates. This program uses a plane wave basis set to perform calculations on periodic systems. All the calculations were performed within the GGA approximation using the PBE functional [35] and adding dispersion corrections [36] (employing the DFT-D method [37]). Pseudopotentials were generated 'on the fly'. The plane-wave energy cut-off was set at 800 eV; a Monkhorst-Pack grid was used for k-point sampling (8 × 8 × 8 grid, spacing 0.001 Å⁻¹). Geometry optimisation was performed using the BFGS algorithm [38]. Since the force-field is needed for aqueous systems, we require a model for the interactions between water and the alkali earth cations and the nitrate anion. We use the SPC/Fw model for water [39] since it is cheap on computer resources and performs well in conjunction with the carbonate and sulphate anions. For completeness, the parameters for the interactions between alkali metal ions and water are listed in table 3 and those for the water model itself are listed in table 4.

Table 4. SPE/Fw parameters for the water molecule (q_H = 0.41∣e∣, q_O = −0.82∣e∣) Wu et al [39] The water-water interaction is given in table 3.

Interaction type	Force constants		Distances		Angles
${V}_{HO}\left(r\right)\,=\,\displaystyle \frac{1}{2}\,{k}_{2}{(r-{r}_{0})}^{2}$	k2 (eVÅ−2)	45.93	${r}_{0}$(Å)	1.012
${V}_{HOO}(\theta )\,=\,\displaystyle \frac{1}{2}\,{k}_{\theta }{(\theta -{\theta }_{0})}^{2}$	kθ (eVrad−2)	3.291			θ0 (deg)	113.24

Molecular dynamics simulations were performed to investigate the high-temperature behaviour of the force-fields. The DL_POLY [40] code was used to perform canonical ensemble calculations over a temperature interval from 1 to 500 K to model lattice expansion. A Nosé-Hoover thermostat with a relaxation time of 0.1 ps was used. The time step was set to 0.5 fs. The selected number of time steps was 50000 (25 ps), of which 10000 (5 ps) were used as an equilibration period. A temperature-scaling interval of 10 K was used, whilst the radial distribution function (RDF) was collected every 1000 time steps, using a bin width of 0.1 Å to obtain the plot. The Ewald summation real space cut-off was set to 8.0 Å, whilst the width of the Verlet shell was 1 Å. The Ewald sum precision was set to a relative error of 10^–5. In order to calculate the long ranged electrostatic (coulombic) potentials the Smoothed Particle Mesh Ewald (SPME [41]) summation method was used. The grid for the k-vector summation was set with dimensions of 8 × 8 × 8, and α, the Ewald splitting parameter, set to 0.36037 Å⁻¹. Isothermal-isobaric (NPT) ensemble calculations were performed using the Nose-Hoover NPT ensemble, with the pressure set to 1 atm and the temperature changed from 300 K to 500 K. The thermostat and barostat relaxation times were 0.1 ps. Other parameters were set as for the NVT simulation.

To obtain the free energies of ion solvation we performed thermodynamic integration (TI) calculations with the LAMMPS code [42]. To perform a thermodynamic integration we must define a pathway between the two systems of interest, in this case the ion in solution and the ion in vacuum. The pathway must be continuous and reversible. The point along the pathway may be defined by the parameter λ. The free energy may then be obtained by integrating over this pathway.

$$\begin{eqnarray}&&{\rm{\Delta }}A={\displaystyle {\int }_{\lambda =0}^{\lambda =1}\left\langle \displaystyle \frac{\partial U\left(\lambda \right)}{\partial \left(\lambda \right)}\right\rangle }_{\lambda }d\lambda \end{eqnarray} \tag{ 3 }$$

where U(λ) is the potential energy of the system at the point along the pathway given by λ. Hence the term in the angle brackets represents the average of the derivative of the potential energy with respect to λ over the course of a simulation at that λ value. We do not use λ directly but instead use a nonlinear function of λ

$$\begin{eqnarray}&&f\left(\lambda \right)=\,{\lambda }^{5}\left(70{\lambda }^{4}-315{\lambda }^{3}+540{\lambda }^{2}-420\lambda +126\right)\end{eqnarray} \tag{ 4 }$$

The use of the nonlinear function stabilises the calculation of the system as λ approaches 0 or 1 [43]. Using the f(λ) of (4) is equivalent to sampling more heavily near the tail end values of λ, with a reweighting scheme given by the derivative of f(λ) with respect to λ. The derivative in equation (3) is obtained via the chain rule

$$\begin{eqnarray}&&\displaystyle \frac{\partial U\left(\lambda \right)}{\partial \left(\lambda \right)}=\displaystyle \frac{\partial U\left(f\left(\lambda \right)\right)}{\partial f\left(\lambda \right)}\displaystyle \frac{\partial f\left(\lambda \right)}{\partial \left(\lambda \right)}\end{eqnarray} \tag{ 5 }$$

where the derivative of the nonlinear function with respect to λ is obtained analytically. The derivative of the potential energy with respect to the nonlinear function is obtained numerically by re-analysing the simulation trajectory and recalculating the potential energy with slightly perturbed values of f(λ) and entering the result into a central differencing scheme for the computation of the derivative. The perturbed values of f(λ) are chosen to be 1% of the value of f(λ) to avoid unphysical negative force-field parameters arising in the central differencing scheme.

The pathway we choose between the solvated and desolvated ion is to deactivate all inter-molecular interactions between the ion and the water simultaneously. We construct a cubic simulation cell containing 3,200 water molecules and equilibrate in an NPT ensemble at 300 K and 0 atm pressure for 1 ns with lattice vectors averaged every 100 fs to get average lattice vectors for the correct density of pure water. We then insert a single ion into the simulation cell (either K or NO₃). This leads to a small additional pressure in the cell when the interactions of the water with the ion are at full strength but is negligible at the simulation cell sizes we use.

The TI simulations were each run in an NVT ensemble for 1 ns with the trajectory printed for re-analysis every 100 fs. We ran 42 simulations at a λ interval of 0.025 (except below the point at 0.025 where the simulations became unstable, however this represents a small error in our calculation and may be excluded). The integration was performed using the trapezoidal rule to obtain the free energy of solvation.

Comparison of our results with 'experimental' tabulations of free energies of solvation is not a simple matter. In addition to ensuring that the results are converged with respect to the size of the simulation cell, it is also necessary to correct for polarisation effects and for the electrostatic potential generated at the site of the ion in a manner consistent with the assumptions that underly the tabulations we wish to compare with. Such issues have been discussed by various authors [44–47] who have reached apparently different conclusions. We consider this matter in the next section.

To demonstrate the reliability of the force field, we have calculated the free energy of hydration of the potassium and nitrate ions using the thermodynamic integration method described in the Methodology section. We must first consider the data which we are to use for the comparison. It is not possible to use experiments to measure the free energy of hydration of isolated ions; what can be measured is the thermodynamics of neutral combinations of ions. This can be used to infer the hydration free energies of isolated ions provided that the hydration free energy of one species can be fixed. The species chosen is usually the proton. So-called 'conventional' free energies are obtained by setting this quantity to zero. However, the free energies obtained by this route cannot be directly compared with simulations because of the arbitrary nature of this assumption. Hence a physically-based model must be used, usually based on the Born model of hydration. These issues are discussed in depth from the experimental perspective by Marcus [46, 47, 49], who is the author of the most extensive tables of thermodynamic parameters for the hydration of ions and more recently by Hofer and Hünenberger [50] from a theoretical perspective. Marcus [47] fixed the absolute free energy of a hydration for a proton at 298.15 K as −1056 ± 6 kJ mol⁻¹. It is now generally agreed from both theoretical [50] and experimental [51] research that a more accurate value is −1100 ± 14 kJ mol⁻¹. Marcus' tabulated values must therefore be corrected by −44z kJ mol⁻¹ where z is the charge of the ion. So his free energy of hydration for potassium (−295 kJ mol⁻¹) should be corrected to −339 kJ mol⁻¹ and his free energy of hydration for the nitrate ion (−300 kJ mol⁻¹) should be corrected to −256 kJ mol⁻¹.

Corrections must also be made to the simulation result. This is true whether spherical or periodic boundary conditions are used [52]. For periodic boundary conditions, the simulation of free energies of ion hydration uses a charged simulation cell, producing a divergence in the electrostatic summation which must be removed. The simplest way of doing this is to add a uniform neutralising correction to cancel the charge. This involves two types of correction. The first, (due to Fuchs [53]), corrects for the interaction of the charges with the neutralising background. The second (first discussed in [54], but usually called the Wigner correction) corrects for the self-energy of the added charge with its periodic images. Most codes using P³M (particle-particle-particle-mesh) approaches (such as LAMMPS, AMBER, DL_POLY) already include these terms, although this may not be immediately apparent (see [55] especially appendix A). Such codes do not include corrections for the artificial polarisation effects caused by the periodic boundary conditions. These are discussed in [56, 57]. Moreover, to connect with the absolute values for individual ions obtained by Marcus and others, we must ignore any effects due to the transfer of the ion through the interface between the gas and water phases and subtract any effects due to the errors in the electrostatic potential at the site of the ion generated by the summation convention used for the electrostatics (see [51, 57, 58] for a full discussion and the importance of the distinction between 'real' and 'intrinsic' quantities when comparing with Marcus' tables). There are also some small corrections arising from choice of thermodynamic conventions [47, 57]. Details are given in the Supplementary Information.

The raw thermodynamic integration results give −297 kJ mol⁻¹ for the free energy of hydration of the potassium ion (the same as the value obtained by Joung and Cheatham [26] using the SPC/E model for water^♯). A hydration free energy of −311 kJ mol⁻¹ was obtained for the nitrate ion. When the corrections are applied, these values become −372.5 kJ mol⁻¹ (corrected Marcus value −339 kJ mol⁻¹) for the potassium ion and −217.5 kJ mol⁻¹ (corrected Marcus value −256 kJ mol⁻¹) for the nitrate ion. However, it is important to remember that experimentally accessible quantities (such as the free energy of solution) depend only on the sum of the free energies of hydration of individual ions. Both the simulation corrections (arising from the electrostatic potential for the cation and anion) and the experimental corrections (arising from the choice of the free energy of hydration of the proton) cancel out. The simulated value for the sum of the free energies of cation and anion is −585 kJ mol⁻¹. The tabulated value obtained from Marcus is −595kJ mol⁻¹ - within about 2% showing excellent agreement.

Further tests using the force-field in water have been performed. Starting from a box of water and potassium nitrate, molecular dynamics simulations using DL_POLY in an NVT ensemble have been performed (2000 water molecules and 50 potassium nitrate formula units). The radial distribution function (RDF) for KNO₃-water was extracted and analysed (see figure 3). The first maximum in the RDF for the nitrogen in the nitrate ion and the oxygen in the water molecule (r_max (N–O_w)) is at 3.59 Å, which shows good agreement with those found in literature, as can be seen from table 10.

The number of water oxygen atoms around the nitrate nitrogen has also been calculated, by integrating the RDF. The coordination number obtained (9.6) shows a good agreement with the values from previous literature as it is possible to see in table 11.

We finally turn to a more detailed consideration of the structures predicted by the simulations. NPT molecular dynamics simulations were performed as previously described and showed that for the lithium nitrate system the structure is calcite-like (space group R $\bar{3}\,$c) at all temperatures and peaks retain their shape (figure 4). The structure loses some long-range order when the system is heated up from 100 K to 500 K (peaks between 4.5 and 5 Å and between 5 and 7.8 Å merge). However, changes in the polymorph are not visible (lithium nitrate is monotropic as experimentally reported). If the temperature is raised to 600 K the long-range order decays further (peaks shift to the right/broaden and their height decreases). At greater atomic separations, the lack of coherence becomes more evident with more peaks merging (around 6, 7 and 9 Å). The behaviour is similar for sodium nitrate (figure 5). At all temperatures the structure is calcite-like (space group R $\bar{3}\,$c); the peaks lose their shape as the temperature rises.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Potassium nitrate is a more interesting case to study, because the RDFs show evidence for three different polymorphs (see figure 6): the ferro-electric trigonal structure (space group R3m) at 100 K; the aragonite structure at 300 K and the trigonal R $\bar{3}$m structure, (quasi-calcite form) at 500 K. On cooling down the system returns to the original structure but not directly. Cooling down the 500 K trigonal structure gives a distorted aragonite structure, whereas the 300 K structure returns to the ferro-electric trigonal form on further cooling. This suggests that the force field is able to reproduce the polymorph changes and shows that those changes are temperature-dependent (potassium nitrate is enantiotropic). Figure 7 shows the good agreement between the experimental K-O distances obtained by X-ray diffraction and those obtained by simulation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size