Understanding the Complex Surface Interplay for Extraction: A Molecular Dynamics Study

Abstract By means of classical molecular dynamics simulation the interfacial properties of methanol and n‐dodecane, which are two potential candidate solvents for use in non‐aqueous liquid–liquid extraction, were assessed. The question of how the interface changes depending on the concentration of extractant (tri‐n‐butyl phosphate) and salt (LiCl) is addressed. Two different models to represent systems were used to evaluate how LiCl and tri‐n‐butyl phosphate affect mutual miscibility, and how the last‐named behaves depending on the chemical environment. Tri‐n‐butyl phosphate increases the mutual solubility of the solvents, whereas LiCl counteracts it. The extractant was found to be mostly adsorbed on the interface between the solvents, and therefore the structural features of the adsorption were investigated. Adsorption of tri‐n‐butyl phosphate changes depending on its concentration and the presence of LiCl. It exhibits a preferential orientation in which the butyl chains point at the n‐dodecane phase and the phosphate group points at the methanol phase. For high concentrations of tri‐n‐butyl phosphate, its molecular orientation is preserved by diffusion of the excess molecules into both the methanol and n‐dodecane phases. However, LiCl hinders the diffusion into the methanol phase, and thus increases the concentration of tri‐n‐butyl phosphate at the interface and forces a rearrangement with subsequent loss of orientation.


Introduction
Hydrometallurgy involving aqueous chemistry is today the most common approach for the recovery of metals,a nd thoughi tw ill remain af undamental toolf or the extraction of many elements, new routes should be explored.I nf act, hydrometallurgy usually involvesacombinationw ith pyrometallurgy, [1][2][3] which was provent ob ei nsufficient for treatingl owgrade ores or residues in an economic way and are weakly selective. The aforementionedr easons, together with the objective to establish ac ircular economy, [4] led to the development of the innovative concepto fs olvometallurgy. [5] The idea of solvometallurgy is to exploit processes similart o those of hydrometallurgyb ut withouta na queous phase. Here, the term "nonaqueous solvent"i su sed in an inappropriate manner,s ince it does not necessarily imply anhydrous conditions, but rather as olvent in which the water content is lower than 50 vol %. This opensabroad range of solvents to choose from,i ncluding molecular organic solvents, ionic liquids, deepeutecticsolvents, and inorganic solvents. [6][7][8][9][10][11][12] Some examples of solvometallurgical processes are the recovery of copper from chrysocolla, [13,14] rare earths and other metals from complexs ilica-rich ores, [15] uranium from carbonate ores, and reprocessingo fs pent nuclear fuel. [16] As in hydrometallurgy, alsoi ns olvometallurgy conventional solvente xtraction is used. In this case though,t he aqueous phase is replacedb yanonaqueoussolvent. [17,18] Although thesep rocesses can be carriedo ut easily on al aboratory scale, on an industrial scale they may be challenging and several conditions should be fulfilled. [5,19] Amongt hose,i n solvent extraction the transition region between the immiscible liquid phases is of paramount importance,a si tc an either facilitateo rh indert he migration of the target compound betweenthe phases and even increase or decreasethe selectivity. For these reasons, overthe years several studies have been undertaken to characterize the interface between solvents for classical solvente xtraction. For example Wipff et al. studied and characterized the interface between aqueous phases and different organic compounds, also in the presence of acids or extractants. [20][21][22][23] They proposed an extractionm echanism that involvest he adsorption of the liganda tt he interface followed by as eries of complexation equilibria between ligandsa nd extraction targets resulting in desorption of the adsorbed complex into the organic phase.T hey also pointed out that desorptiono ft he complex might be facilitatedby an increased interfacial concentration of the complex, the extractant molecules,o ra ny other surfactant, whichw ould reduce the interfacial pressure. [21] Moreover,s ince the complexationo ft he target by the ligand seems to occur at the interface, the study of ligand affinity towards pecific targets( e.g.,b yq uantum chemicalm ethods) should not neglect the impact of the interface structure, and therefore ad eeper knowledgeo ft he structural features of the interface is mandatory.
Most of the previous studies were focused on providing fundamental insight into designing aqueous solvent extraction. Obviously,i nv iew of the emerging concept of solvometallurgy, these principles should be expanded to ab roader set of immiscible pairs of solvents. Previously,w einvestigated the mixture of methanol( MeOH)a nd n-dodecane (DD) in order to discuss phase separation, and studied how the system can be altered and the leadingp rinciples for optimizing the system. [19] In the present work, we employed the same mixture to characterize the interface between two immiscible molecular solvents by means of molecular dynamics simulations.W ea ssess at a molecular level to what extent the solvents form distinct phases and study the structuralf eatures of the interface between them.I no rder to do so, we also added as urface-active ligand,n amely tri-n-butylp hosphate (TBP), and as alt, namely LiCl, which was shown to improvet he phase separation of these solvents, to the system.T his is af irst step towards the understanding of how the interface structure can influence solvent extractionp rocesses and therefore of how thesep rocesses can be improved.

Investigated systems
We simulated 15 systems, which can be divided into two categories. The first category contains eight biphasic systems simulated in cubic boxes (Figures 1a nd 2). These systems were fundamental first to assess how well the force field could represent the phase separation, and second to analyze the mutual miscibility of the solvents at different ligand and salt concentrations. Therefore, their initial configurations were randomly generated to ensure that the obtained phase separation and the interface between the formed phases were not conditioned by the starting geometry or other constraints. The second category contains seven systems, for which the simulation box is longer in the z direction ( Figure 3). These simulations allow us to investigate interfacial structure and the behavior of the components adsorbed on the interface (technically,t wo interfaces are generated because of periodic boundaries conditions, but in order to have control over the analyses, TBP was placed only at one interface). Different from the first category,i n which the cubic simulation boxes were packed with arandom configuration, systems belonging to the second category were generated by juxtaposing the two phases, in order to facilitate equilibration. Thus, the generated interfaces are flat and the behavior of the adsorbed molecules can be easily studied. The first and second categories of systems are labeled with the initial strings CUB-and NCUB-respectively.I nb oth cases the remaining letters of the labels represent the components of the system, that is, Ms tands for MeOH, Lf or LiCl, Df or DD, and Tf or TBP. More details of the compositions of the systems are listed in Ta ble 1.

Computational details
The initial configurations of the systems were generated by using the PACKMOL package (version 17.039). [24] An initial cell vector of  6.4-6.7 nm was used for the CUB systems, depending on the number of molecules in the cell. Molecules were randomly placed within these boundaries. For the NCUB systems, MeOH and DD were placed in equally sized opposing cells of 6.0 nm length in x and y directions and 4.0-4.3 nm length in the z direction. TBP was randomly placed within 0.5-1.1 nm distance to the interface between the two phases, depending on the number of TBP molecules. Li + and Cl À were placed randomly and independent from each other within the MeOH phase. Classical molecular dynamics was performed by using the LAMMPS program package (version 17 Nov 2016). [25] The well-known OPLS-AA force field was used for MeOH, DD, and LiCl, [26,27] whereas for TBP we opted for the force field recently developed by Ali et al. [28] since it was shown to perform very well when TBP is mixed with n-dodecane. Nonbonded interactions were described by the 6-12 Lennard-Jones potential and Lorentz-Berthelot mixing rules for pairs of different atoms. [29] Ac utoffo f1 .6 nm was selected for the calculation of Lennard-Jones and Coulombic interactions together with ap article-particle particle-mesh solver mapping the atom charge to a3 Dm esh. [30] Equilibration of the systems was obtained by simulating for 2.3 ns with NVE, NpT,and NVT ensembles. After an initial energy minimization, the systems were simulated for 0.3 ns in an NVE ensemble with added velocity scaling, during which the simulation box was deformed so as to reach ap re-equilibration density of 0.8 gcm À3 . Following this, the systems were simulated for 0.5 ns in an NpT ensemble by applying the NosØ-Hoover chain thermostats and barostat to achieve constant pressure and temperature (T = 297.15 K, t = 100 fs and p = 1bar, t = 1000 fs, respectively). [31,32] Finally,u sing the average volume of the previous 0.05 ns as final box volume, the system was equilibrated for another 1.5 ns in an NVT ensemble with the same settings for the thermostat. The subsequent production run consisted of 10 ns of simulation time in an NVT ensemble, with the same settings as during equilibration. The time step was set to 0.5 fs during the whole procedure, and every 1000th time step was saved in atrajectory for further processing.

Computational analyses
The obtained trajectories were analyzed with the TRAVIS software and other in-house scripts. [33,34] TRAVIS offers different kinds of functions allowing the analysis of the interaction among the components of the systems. Intra-and intermolecular interactions can be taken into account.
The domain analysis feature implemented in TRAVIS, which is based on Voronoi tessellation, was employed to study the phase separation of the solvents and mutual miscibility. [35] Radical Voronoi tessellation [35][36][37] was used on every saved trajectory time step to calculate the number of molecules that migrated from one phase to the other,a nd thereby to study how mutual miscibility changes in different conditions. The Voronoi tessellation works as follows: the system is divided into subsets, wherein as ubset is as et of specific atoms. Subsets can either contain all atoms belonging to a molecule or atoms selected from different molecules depending on the aim of the analysis. Once the atoms are flagged for the subset they belong to, the algorithm analyzes the interconnection of the atoms belonging to the same subsets. Atoms of the same  subset that are in close contact with each other are considered part of the same domain, and the total number of domains formed by each subset describes how the subset is distributed in the system. When the subset contains all atoms of as pecific compound, an umber of domains equal to unity means that this compound tends to aggregate. On the other hand, an umber of domains equal to the number of molecules of the compound means that the compound is fully solvated by another solvent. Using radical Voronoi tessellation, Brehm and Sebastiani were able to successfully determine the exposed surface of 1-ethyl-3-methylimidazolium acetate droplets in vacuum and used it to describe a liquid-vacuum interface. [38] In this work, which treats al iquidliquid interface, raw data from the Voronoi tessellation analysis were processed to obtain information on the number of molecules of one solvent or compound dissolved into the other.T his was possible because of the neat phase separation of the two solvents. Phase separation generates two main domains that contain most of the atoms belonging to the subset. Therefore, by counting the number of atoms of the subset that do not belong to the main domain, it is possible to calculate the number of molecules that migrated into the other solvent. For the Voronoi tessellation-based analyses in this work, MeOH and DD were treated as separate subsets, while TBP and LiCl were included as Voronoi sites but were not analyzed further.
Combined distribution functions were generated by correlating the positions of TBP molecules with their orientation and by applying the Kernel Density Estimation provided in the Python library Seaborn. [39] The orientation of the molecules was evaluated by the angle a between the z axis and the vector whose base and tip are positioned on Pa nd Oo ft he P=Ob ond, respectively,a ss hown in Figure 4. The 1/sin a cone correction to correct the uniform angular distribution was applied. [33] Results

Mutual dissolution of solvents
The distribution of the componentsi nt he different formed phases was analyzed with the CUB systems. We used Voronoi tessellation (see Methods Section) to calculate the number of molecules that migrated from one phase to the other in each time step. The resultsa re summarized as histogramsi n Figure 5, which show how often acertain number of molecules were dissolved in the opposite phase within the simulation time. The distributions depict very clear trends depending on the selected conditions. The left panelso fF igure 5s how the histograms related to the number of MeOH molecules dissolved in DD. The top and bottomp anels refer to the cases withouta nd with LiCl, respectively,w hile colors refer to different concentration of TBP in the systems. To evaluatet he influence of LiCl on the mixing, we comparethe width of the distributionsa nd the position of the maxima. In the presence of LiCl, the number of molecules in the opposite phase is decreased. This is especially clear when comparing the right tails of the distributions:i nt he LiCl-free cases the tails extends to high numbers, up to 24 (not shown in the figure) at aTBP concentration of 60 molecules, whereas it stops at ten MeOH molecules in the presenceo fL iCl. Since in industrial extraction processes neat phase separation with low mutual miscibility of the solvents is mandatoryt oa void solvent losses, according to these analyses the addition of LiCl might be as olution in cases in which the mutual miscibility of the solvents is too high. The same evaluation was done for different TBP concentrations. The neat systems (red color) exhibit the lowest dissolution of MeOH molecules into DD, whereas the increased TBP concentration leads to enhanced migration of MeOH into DD. The maximum number of migrated MeOH molecules is found at the highest TBP concentrationi nt he absence of LiCl. As mentioned above, mutual miscibility of the solvents should be avoided; therefore, according to the analysis, for this specific solventp air the amounto fT BP used for the extraction should be carefully evaluated to minimize solventlosses.
The right side of Figure 5s hows the same kind of analysis as the left butw ith respect to DD molecules migrated into MeOH.I ng eneral, the number of dissolved DD molecules is lower than the number of MeOH molecules dissolved in DD. The comparison of the top-right and bottom-right panels shows the strong effect of LiCl. Its presence hinders the migration of DD molecules into MeOH,i nl ine with what we found in our previous work. [19] The effect of TBP was less strong;n evertheless, we point out that, as expected, for higher concentrations of TBP the migration seems to be facilitated, as shown by the barso nt he right sideo fb oth panels. Hence, with regardt ot he migration of DD into MeOH,L iCl might againb e used to decreasem utual miscibility,w hereas the amount of TBP should be carefully evaluated.
While the data depicted in Figure 5i llustrate distinct trends, it might be interesting to study the temporal evolution of the number of solvated molecules. These plots may be found in the Supporting Information.
To further analyze the migration of molecules into the other solvent, Ta ble 2l ists the average numbers n MeOH and n DD of MeOH molecules solvated in DD and vice versa. In the case of MeOH ani ncrease in n MeOH can be observed on addition of TBP.A lthough no significant difference in migration is visible for systems with 15 or 30 TBP molecules, increasing the TBP concentration to 60 molecules leads to an average number of solvated molecules twice as high. Interestingly,t he TBP concentration has no significant effecto nD Dm igration. However, the addition of LiCl leads to as ignificant decrease in migration for both MeOH and DD, in accordance with our earliero bservations. We note that, since the domain analysisincluded TBP as Voronoi sites, it is possible that molecules solvated by TBP rather than the opposite solvent are counted toward the number of migrated molecules. However,i fT BP is excluded, molecules solvated by the opposite solvent and separated by TBP might be counted toward the main domain. It is cleart hat the true number of migrated molecules must lie in between. Hence, we carried out additional analyses that excluded TBP as Voronoi sites, which can be found in the Supporting Information. Althought he number of migrated molecules decreases, similar trendsa re observed, that is, the observed effects on mutual miscibility are not merely due to solvation by TBP.
Next, we turn to the systemst hat were set up with an interface (see Figure3). Figure 6s hows the normalized density profiles of the solvents, centered at the interface calculated for NCUB systems. An ideal biphasic system would show an eat and sudden drop of the density of the solvents at the interface position. Changes in the concentrationo fT BP,i ndicated by differentc olors in the plots, strongly affect the behavior of the solvents at thei nterface. We focus first on the curvesd escribing MeOH density (the curvess tarting at about 2o nt he right), as the noise in the DD curvesm akes their evaluation complicated.F or higherc oncentrations of TBP the density curves of each solvento nt he solventp hase side were lower,w hereas they were higher on the other side. As observed in the previous analyses, this effect was counteracted by the addition of LiCl, as it decreases the densities of the solvents in the opposite phase and increases the densities of the solvents in their own phase. This means that TBP increases mutual miscibility in the vicinity of the interface, whereas LiCl can be used to reducem utual miscibility.T his interfacial exchange is in line with the mixing behavior observed in the Voronoi analyses and it provides ag ood alternative toolt ostudy solventm iscibility from ad ifferenta nd complementary perspective, as it provides information on the behavior of the solvents at the interface.

Interfacestructure
In the previous section it was determined how TBP concentration affects the mutual miscibility of the solvents by increasing the exchange of molecules between the two phases.T herefore, the positioning of TBP molecules was studied by visual analysisoft he trajectories of the CUB systems. In Figures 1a nd  2t wo snapshots of the systemsa fter equilibration are depicted, which represent the general situations over the whole simulation time. Distinctively biphasic systemsw ere obtained. From visual inspection it is apparent that in both cases the TBP moleculesa re adsorbed at the interface. Therefore, af urther investigation of TBPa dsorption at the interface seems of paramount importance. With this aim, we turn again to the NCUB systems. Figure 3s hows the simulation box of system NCUB-MDT30. As expected TBP was mainly adsorbed at the interface betweent he solvents, and phase separation was maintained during the whole simulationt ime. From the snapshot we also obtained afirst glimpse of the possibility of desorption of TBP into each solvent. Figure 7s hows the density profiles of solvents and TBP along the z direction. The solid red curves represent TBP density alongt he z direction in systemsw ithout LiCl, and the dashedr ed lines the TBP density in presence of the salt. The green solid line represents the density along the z direction of LiCl, and black and blue lines describe the density along the z directiono fM eOH andD D( i.e.,t he MeOH phase starts from 2000 pm, reaching into the DD phase, and vice versa),r espectively.The three panelsr eport information for different concentrationso fT BP,w hich was increased from top to bottom. TBP exhibits very different behaviors depending on its concentration and the presence of LiCl. For increasing concentrations of TBP,w hich was shown before to result in an egative effect on phase separation, we find broader density profiles, suggesting the mixingo fT BP into MeOH and DD, at least in the vicinity of the interface. LiCl, on the other hand, tends to increase the density of TBP in the interfacialr egion. This can be explained by analyzing the shapeo ft he density curvesi nt he surrounding of the interface. We observe that LiCl decreases the density of TBP on the MeOH side, whereas it increases on the interface and remains the same on the DD side. Hence, this effect could be relatedt oah indrance of the diffusion of TBP molecules into MeOH due to LiCl, the subsequent accumulation of TBP at the interface, and, when the interface is saturated, into DD.
To further study the structural features of the interface, we show in Figure 8t he correlation between the positioning and the orientation of TBP molecules. The x axis represents the position of TBP molecules during the simulation, along the z direction definedb yt he position of the Pa tom, while the orientation, which is defined as the angle a betweent he z axis of the simulation and the P=Ob ondo fT BP (see Figure 4) is reported on the y axis. Further visualizations of this analysisi ncluding cosa and the order parameter S ¼ 3cos 2 aÀ1 2 can be found in the Supporting Information. Additionally,s imilar analyses of the MeOH orientation in dependence on its position along the z axis were carriedo ut and can be found in the Supporting Information with as et of combined distribution functions of the hydrogen-bond angle and distance between MeOH and TBP.F or values of a close to 1808,t he orientation of the molecule has the polar moiety pointingt oward the DD  phase, and for values close to 08 an orientation toward the MeOH phase is indicated.A se xpected, TBP is found preferentially at the interfacial region between the two liquids (at approximately 60 000 pm), whichc onfirms the above observation that it is adsorbed at the interface. Figure 8a lso allows the comparison of interfacialT BP orientation depending on its concentration and the presence of LiCl. In the case without LiCl, TBP occupies aw ider portion of space,a nd diffusioni nto both phases seems possible. When TBP is adsorbed at the interface, it exhibits ap referentialo rientation with a 408 in which the polar group points at the MeOH phase and the butyl groups at the DD phase. On the DD side of the interface, it is apparent that TBP lacks preferentialo rientation. On the MeOH side, however, orientations closer to 08 still seem to be slightly favored. As shown by the panels on the right side of Figure 8, the situation radically changeso nce LiCli sa dded. In fact, in the presence of LiCl, the TBP molecules are constrained to the interface on the MeOH side, whereas they still exhibit a certain degree of freedom on the DD side. Even more interesting is the effect of this constraint on the orientation of the molecules at the interface. We observe that, on increasing the number of TBP molecules, a is more evenly distributed between 0a nd 1808,p robablyd ue to an excesso fT BP molecules at the interface, which can be solvedo nly by rearranging the interfacial structure itself. The mostc ommon orientation in the presence of 60 TBP molecules of a % 808 depictsacompletely different situation, in which TBP is aligned parallel to the interface. The most reasonable interpretationo ft hese analyses might be that, with the addition of LiClt ot he polar phase,t he excesso fT BP molecules cannotd iffuse into the MeOH phase. This is also evident in the reduced interaction between MeOH and TBP after addition of LiCl, which is visible in the combined distribution functions reportedi nt he Supporting Information. Hence, TBP loses its preferentialo rientation towards the solvents, so that its P=Obond aligns with the xy plane and the interactions with other TBPm olecules are maximized. This is important with regard to solvent extraction, because the lack of specific TBPo rientation in systemsc ontaining LiCl most likely negativelya ffects al igand-target extraction process. In fact, TBP would not expose the active side (the phosphate group) to the polar solvent, which is where the target is usually dissolved. On the other hand, for low concentrationso fT BP the orientation seems to be maintained also in presence of LiCl, and therefore ab alance betweent he concentration of LiCl and TBP might be the key to obtaining good phase separation withoutimpeding the extraction process.

Conclusion
We have reportedm olecular dynamics simulations of the interface between MeOH and DD in the presence of LiCl and TBP, which is ac ase study fors olvent extraction in solvometallurgy. In the current work, we consideredo nly one concentration of LiCl. Different salt concentrationsm ight be the subjecto ff urther studies, since the concentration of LiCl must be carefully balanceda gainst potentially detrimental effects such as reduced metal solubility in the polar phase or decreasing extraction efficiencyb yc ompetition of Li + with the extraction process of another metal ion. Additionally,s alts other than LiCl may shows imilare ffects andm ight be worthy of consideration. Indeed, our previous work suggests that the increased phase separation is not directly related to the nature of LiCl, but rather stems from as trengthening of the hydrogen-bond network in the methanolp hase. [19] Thus, we expect similar results foro ther alkali metal halides. For industrial application of solvente xtraction,o ne of the main requirements is neat phase separation between the solvents. Thus, with the cubic model of these systems and an ewly developed analysis based on the Voronoi tessellationm ethod, we evaluated the effect of TBP and LiCl on the mutualm iscibility of the solvents. We found that higher TBP concentrations lead to enhancedm utual miscibility,a nd that LiCl counteracts this effect and improves phase separation. We observed that TBP is mainly adsorbed at the interface betweent he two solvents. Therefore, we focused on the interface andi ts features related to TBP and LiCl. With this aim, we built seven noncubic systems, startingw ith af lat and well-fixed interface between the solvents in as pecific position of the box. These simulations allowed us to study how the structuralp roperties of the interface (relatedt oa dsorbed TBP molecules) vary depending on the chemical environment. By plottingt he correlation between the positiono fT BP along the z axis and its orientation,w es howed that TBP adsorption changed depending on its concentrationa nd LiCl presence. Without LiCl, adsorbed TBP exhibits aw ell-definedo rientation, with its polar moiety oriented towards the MeOH phase and its alkyl chains pointingtowards the DD phase, whichmaximizes the polar-polar interactions on one side and the nonpolarnonpolari nteractions on the other.F or higherT BP concentrations this feature was preserved, and the TBP molecules that could not be adsorbed in the first layer at the interface owing to lack of space diffused into the solvents.
The presence of LiCl, which slows down the dynamics of the MeOH phase (as described in our previous work), [19] hindered the diffusion of TBP molecules into the MeOH phase and increasedt he concentration of TBP at the interface. Also the orientation of interfacialT BP was affected, especially for higher concentrations of TBP.I nf act, with LiCl, when the TBP concentration was increased, TBP migration into MeOH was hindered, and this constraint led to ar earrangemento ft he TBP layer,i n which TBP molecules oriented towards each other to maximize the interaction of the polar and nonpolar moieties. All these effects can affect solvent extraction. In fact, in the absence of LiCl, the increased exchange of interfacial TBP molecules and a well-defined orientation of TBP molecules most likely positively affect the formation of al igand-target complex and its migration from one phase to the other.O nt he other hand, TBP was provent oi ncreaset he mutual miscibility of the solvents, which should be avoided in industrial application.I ndeed, the works of Wipff et al. suggest that asaturated interface is necessary for extraction to occur,w hich limits the range in which TBP concentration can be optimized against detrimental effects such as solvent mixing. The addition of LiCl can improve the system by decreasing mutual miscibility,b ut it also interferes with the orientation of TBP,w hich might affect complex formation,and therebyh inderthe extraction.