Molecular Dynamics Simulation of Wetting and ... - ACS Publications

Article Cite This: J. Phys. Chem. C 2019, 123, 13551−13560

pubs.acs.org/JPCC

Molecular Dynamics Simulation of Wetting and Interfacial Properties of Multicationic Ionic Liquid Nanodroplets on Boron Nitride Monolayers: A Comparative Approach Farhad Ghalami,† Tahereh Sedghamiz,‡ Elaheh Sedghamiz,*,† Fatemeh Khashei,† and Ehsan Zahedi† †

Chemistry Department, Shahrood Branch, Islamic Azad University, Shahrood, Iran 36199-43189 Medicinal and Natural Products Chemistry Research Center, Shiraz University of Medical Sciences, Shiraz, Iran 713485-3734

‡

Downloaded via BUFFALO STATE on July 31, 2019 at 07:36:01 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: In the present work, wetting and interfacial properties of two new kinds of imidazolium-based ionic liquid (IL) nanodroplets, that is, dicationic IL and tricationic IL with respect to a baseline of a traditional monocationic IL on a boron nitride (BN) surface, is studied using classical molecular dynamics simulation. Calculated contact angles of all studied nanodroplets on the BN surface are lower than 90°, which demonstrate favorable wettability of the studied ILs on the BN sheet. Density profiles show an enhancement in density of both cations and anions at the BN surface for all studied droplets. Orientational analysis revealed that for monocationic ionic liquids, the imidazolium ring of cations lies parallel to the BN surface, while for multicationic ILs (dication and trication), the perfect parallel arrangement could not be observed. Identification of the truly interfacial molecule method has been used to analyze the structure and dynamics of the ions at the interface. It shows that the survival probability of cations at the interface decreases with increasing number of imidazolium moieties in the cation. The above observations signify the importance of charge and size of the cation on wetting and interfacial properties of ILs in contact with solid surfaces.

1. INTRODUCTION Investigation of small volumes of liquids on solid surfaces has attracted considerable attention over the last decade thanks to their potential applications of relevance to energy conversion, surface coatings, sensors, drug delivery, and biological applications.1,2 The shape of the drop resting on the surface depends on the molecular properties of the liquid droplet and the characteristics of the surface.3 The equilibrium balance of the forces acting on the liquid−vapor interface in contact with a solid surface determines the so-called “equilibrium contact angle”, which is a useful measurement to quantify the interfacial energy between solid, liquid, and vapor phases.4,5 Wetting behavior of water droplets on different surfaces has been widely studied using both experimental and theoretical approaches. These studies have shown that the wetting and dynamical properties of water can be strongly dependent on the chemical composition and microstructure of the contacting surfaces.6−8 However, unavoidable impurities and defects of the surface, which essentially changes the liquid/solid interfacial properties, can limit the experimental analysis of these systems.9 For instance, experimental observation indicates that the water contact angle on a graphite surface at room temperature ranges from 60 to 80° because of the sensitivity of contact angle measurement to impurities or defects on the surface.10 Classical MD simulations have been © 2019 American Chemical Society

widely used to study the wetting properties of different solid surfaces since the pioneer work of Rahman and Stillinger.11 Although quantum mechanical calculations are more precise compared to classical molecular dynamics (MD) simulation, they are limited to the system size. A system containing thousands of atoms could be investigated by classical MD simulation using empirical force fields with an acceptable consistency between the predicted properties based on classical MD simulation and experimental data.12,13 Among nonaqueous media, interests in ionic liquids are increasing, especially in surface wetting applications, due to their favorable conductive and vapor pressure properties and most importantly the possibility of designing task-specific ILs through the suitable combination of cations and anions.14 Malali et al. investigated the wetting behavior of an ionic liquid droplet containing 1-methyl-3-butylimidazolium and fluorophosphate on the titanium dioxide (110) surface using molecular dynamics simulations. They found that the imidazolium rings of cations lied parallel to the titanium dioxide surface, but their orientation in the upper layer was almost vertical.15 Castejón et al. implemented a fully Received: December 13, 2018 Revised: May 7, 2019 Published: May 8, 2019 13551

DOI: 10.1021/acs.jpcc.8b11987 J. Phys. Chem. C 2019, 123, 13551−13560

Article

The Journal of Physical Chemistry C

more probable, which is the result of formation of a strong network of hydrogen bonds within a trication. Slab-geometry MD simulations revealed that the behavior of trications at the liquid/vacuum interface is completely different compared to those observed for cations of traditional ILs. Furthermore, orientation of trications at the interface was highly dependent on their alkyl chain length. While LTILs with shorter alkyl chains (i.e., n = 3 and 6) form an inverse-arc shape structure, LTILs with longer alkyl chains (i.e., n = 10) form a sinuous structure at the interface. However, although our previous studies unfolded diverse properties of LTILs in bulk and at the liquid/vacuum interface, efforts are still needed to answer the question of whether the additional cationic moieties in LTILs impose a specific influence on the structural arrangement of multicationic ILs in contact with solid surfaces and therefore exhibit different interfacial properties. Investigating the wetting and interfacial properties of multicationic ionic liquids on a BN monolayer could provide new insight into the behavior of multicationic ILs in contact with solid surfaces and could open the possibility of new applications beyond traditional monovalent ILs. Motivated by these considerations, in this work, we performed classical MD simulations to model monocationic ionic liquid (MIL), dicationic ionic liquid (DIL), and tricationic ionic liquid (TIL) nanodroplets all having [Tf2N]− as the anion (Figure 1) in contact with a monolayer hexagonal BN sheet (Figure 2). In the following, we investigate the effect of multiple cationic groups in imidazolium-based ILs on the contact angle, solid/liquid interfacial structure, orientation, and composition in contact with the h-BN surface for the first time.

phenomenological approach to calculate wetting and transport properties of ionic liquids. They illustrated that, despite the complexity of ionic liquids, they behaved similarly to simple liquids when interacting with a solid surface.3 Wetting behavior of ionic liquid nanodroplets of different sizes on a graphene sheet was studied by Herrera et al. using molecular dynamics simulations. They reported a strongly adsorbed layer on the graphene sheets, which was highly structured, and both ions adopted particular arrangements to maximize the interaction with the graphene. They showed that contact angles increased remarkably from the smallest nanodroplets to the largest ones.16 Due to the remarkable chemical and physical properties of graphene, other 2D nanosheet analogues have also attracted much attention. A graphene structural analogue, hexagonal boron nitride (h-BN), exhibit weaker hydrophobicity due to the polar nature of BN bonds.17 BN nanostructures are also wide band gap semiconductors, good electrical insulators, and highly resistant to oxidation, and therefore they are more suitable for very high temperature (up to 1000 K) applications in which graphene would burn.18 There are few reports that focus on the adsorption properties of ILs on the BN surface. Garcia et al.19 studied the adsorption of choline benzoate ([CH][BE]) ionic liquid on four types of nanosheets, that is, graphene, silicene, germanene, and boron nitride using a density functional theory (DFT) approach. They found that the interaction between the IL and all studied nanosheets mainly featured van der Waals forces, leading to strong benzoate ion−surface π stacking. They also observed that IL− graphene and IL−BN provided similar features such as the interaction mechanism, charge transfer, binding energies, or electronic structure. Shakourian-Fard et al.20 studied the adsorption properties of three types of cations, that is, [Bmim]+ , [Bpy]+ , and [Btma]+ , with different anions ([BF4]−, [PF6]−, [Tf2N]−) on the h-BN surface using quantum chemical calculations. They found that upon adsorption of ionic liquids on the h-BN surface, the overall charge on the cation, anion, and h-BN surface changed and charge transfer between ILs and the h-BN surface occurred. In addition, the contribution of the dispersion component (ΔEdisp) in each complex was more than that of the electrostatic one (ΔEelect). The thermochemical analysis also indicated that the free energy of adsorption (ΔGads) of ILs on the h-BN surface was negative, and thus the adsorption occurred spontaneously. To date, most observations/results in this field are limited to monocationic ILs, that is, with monovalent cations. Symmetric multicationic ILs as a new category of ionic liquids have a cation with two or three identical cationic moieties bridged by alkyl chains. Hence, each cation in multicationic ILs carries two or three positive unit charges, which could achieve higher stability, greater tunability, and lower vapor pressure compared to MILs.21−23 In our previous studies, we investigated the structural, dynamic, and interfacial properties of linear tricationic ionic liquids (LTILs) by means of classical MD simulations and DFT calculations.24−26 MD simulation of LTILs having different chain lengths showed that organization of anions around cations was almost similar to those of MILs and DILs. In addition, the diffusion coefficients of the studied LTILs were smaller than those of both MILs and DILs with comparable viscosities. Also, similar to MILs and DILs, cations had a major role in carrying electric current in LTILs, but this role became bigger from MILs to TILs. DFT calculations revealed that a helical structure for the cations of LTILs is

Figure 1. Molecular structure of the (a) MIL cationic moiety (1propyl-3-methylimidazolium), (b) DIL cationic moiety (1,3-bis(3methylimidazoliume-1-yl) propane), (c) TIL cationic moiety(1-(1′methyl-3′-propylimidazolium)-3-(1″-methyl-3″-propylimidazolium) imidazolium) and (d) [NTF2]− anion. 13552


Article


Figure 2. Snapshot of the model system at initial time of the simulation: (a) IL nanodroplet, (b) BN sheet, and (c) BN−droplet initial configuration. Color code: green, cations; red, anions; blue, nitrogen atoms; pink, boron atoms.

2. METHODOLOGY 2.1. Simulation Details. MD simulations of MIL, DIL, and TIL on the monolayer BN sheet were performed using the DL_POLY 2.17 package,27 and structures were visualized by VMD software.28 Ionic liquid interactions were described using a non-polarizable force field, which is based on the systematic all-atom force fields (AA-FFs) developed by Canongia Lopes and Pádua29 for imidazolium-based ionic liquids with a slight modification by Tsuzuki and co-workers.30 We have used a combination of parameters derived from quantum chemical calculations and the mentioned force field. Atomic charges for ions were calculated using the ChelpG method from previous DFT calculations (reported in Tables S1, S2, and S3) and were used in the force field in the intact form (without charge reduction).30 The accuracy of the used force field for the studied ILs has been established in our previous publications by simulating their dynamic, structural, and interfacial properties.24−26 The parameters of the force field could be found in ref 30. The monolayer hexagonal BN sheet was constructed with 7542 BN atoms with dimensions of 14 × 14 nm, which is kept frozen to prevent it from shifting during the MD simulations. A non-electrostatic model of the BN sheet was employed in all simulations according to the results obtained from previous studies.31,32 Kamath et al.33 employed both electrostatic and non-electrostatic models of the BN surface in the simulations of ILs in contact with the BN surface and reproduced the same overall behavior, observing only minor discrepancies between the computed energetics and structural features for the two different models.34 The van der Waals contributions were included by Lennard-Jones 12−6 (LJ) potential, and the Lorentz−Berthelot mixing rule was used to model interactions between the liquid and the BN sheet. For each system, IL ion pairs (also ion triplets and ion quartets) were first replicated in order to reach a cubic box. In fact, droplet size effects on the contact angle and therefore the number of particles in each droplet were selected carefully in order to finally reach the droplets with the same sizes (Table 1). Afterward, spherical droplets were obtained by performing

Table 1. Specification of the Simulation Boxes for the Studied Droplets on the BN Surface IL

no. of ions (anions + cations)

no. of atoms

droplet diameter (nm)

MIL

432

7992

6.03

DIL

375

7875

6.04

TIL

336

7980

6.10

slab size (nm) x = y = 15, z = 20 x = y = 15, z = 20 x = y = 15, z = 20

heating−cooling (600−298 K) without applying periodic boundary conditions in ∼3 ns for MIL and 5 ns for DIL and TIL. Nanodroplets of MIL, DIL and TIL were placed on the center of a BN nanosheet with a 4 Å initial distance, and periodic boundary conditions were applied in all directions. The size of the slabs was selected to be large enough to allow IL nanodroplet spreading without being influenced by the slab edge planes’ asymmetric configuration. The structures of the BN sheet and droplet and the BN−droplet initial configuration are shown in Figure 2 for MIL. The specifications of the simulation boxes are provided in Table 1. An initial equilibrium was reached allowing nanodroplets to come into contact with the surface by performing a series of MD simulations with the Nosé−Hoover thermostat started at 600 K and subsequently annealed gradually to 298 K in ∼5 ns.35 Afterward, the simulations were continued in 298 K allowing nanodroplets to attain a constant contact angle in 10 ns. Finally, a 10 ns production run was performed under constant NVT conditions at 298 K using the Nosé−Hoover thermostat (relaxation time, 0.2 ps) for the analysis. A shortrange cutoff distance of 5 nm was used for Lennard-Jones and electrostatic calculations. The motion equations were solved by using the Verlet leap-frog integrator with a 1.0 fs integration time step. The Coulombic long-range interactions were calculated using the Ewald method with a precision of 1 × 10−5.36 2.2. ITIM Analysis. To precisely investigate the adsorption behavior of the studied ILs on the BN surface, ITIM analysis can be of great help. The ITIM code has been specially 13553


Article


Figure 3. Illustration of the contact-angle measurement based on MD simulations for (a) MIL, (b) DIL, and (c) TIL nanodroplets on the monolayer BN sheet. The blue dots denote to the liquid/air interface, and the blue curve marks the fitted surface of the nanodroplets. Color code: green, cations; red, anions.

the center of mass of the droplet is defined. The density profiles are obtained for the slabs with a 2 Å width perpendicular to the BN surface. The liquid−vapor interface is defined to be 0.2 g cm−3 at each slab, and thus the radii (r) of the droplet is considered as the distance from the z axis passing through the center of mass and the position of the liquid−vapor interface at each slab. Subsequently, the obtained profile of the droplet radius is fitted to the second-degree polynomial. The contact angle is calculated from the fitting polynomial at a height of 3 Å above the BN surface. Figure 4 illustrates the time evolution of the contact angle for nanodroplets of MIL, DIL, and TIL on the BN surface. It can be observed that MIL and DIL reached the equilibrium contact angle at roughly 4 ns, whereas it is reached at roughly 6 ns for TIL. Therefore, all the contact angle data reported in the following sections were calculated as averages in the 6 to 10 ns simulation time frame. The contact angle error bars for nanodroplets of MIL, DIL, and TIL on the BN surface as a function of time for 10 time points of the last 2 ns of MD simulation are reported in Figure S1. According to Table 2, all droplets adopt contact angles lower than 90°, indicating that liquid−surface interactions are favorable in all of the studied systems. Surprisingly, the contact angle increases slightly from MIL to DIL but increases drastically from DIL to TIL. Furthermore, by increasing the contact angle from MIL to TIL, the height of the nanodroplets increases, while the adsorbed

developed to investigate interfacial properties by providing a precise list of interfacial molecules. The detailed information on methodology of ITIM analysis could be found in the original paper.37,38 In this study, the ITIM code has been used with a slight modification in order to be applicable to use for our model systems. The surface is approximated by the set of intersection points of the probe spheres with the test lines along which they are moved and at the position where the probe spheres were stopped. The distance of two neighboring test lines has been set to 0.5 Å in the yz plane, and the test lines were arranged in 50 × 50 grids along the xy plane of the interface. The radius of the probe sphere has been set to 2 Å for MIL, 3 Å for DIL, and 4 Å for TIL to ensure that it is in the order of the size of the molecules.

3. RESULTS AND DISCUSSION 3.1. Contact Angle. Figure 3 displays the snapshots of the studied droplets on the BN surface at the final equilibrium stage. It shows that all liquids equilibrated to drops with a quasi-hemispherical shape. The degree of wetting of a liquid on a solid surface could be determined by measuring the contact angle of the liquid on the surface. In this study, calculation of the contact angle (θ) has been performed based on the method reported by Giovambattista et al.39 According to the procedure, a z axis, normal to the BN surface, passing through 13554


Article


Figure 5. Number density (ρ) profile for (a) MIL, (b) DIL, and (c) TIL nanodroplets as a function of height on the BN sheet.

Figure 4. Contact angle θ for (a) MIL, (b) DIL, and (c) TIL nanodroplets on the h-BN surface vs molecular dynamics simulation time.

nanodroplets, highlighting lower solidation of the TIL nanodroplet on the BN surface. In all three graphs, the peak of the cation appears slightly ahead of the peak of the anion, indicating that although both ions tend to be adsorbed on the BN substrate, cations have a higher tendency to be adsorbed on the BN surface. Previous MD simulation studies have also shown that cations have a higher tendency to be adsorbed on the graphene surface compared to the anions in imidazoliumbased ionic liquids.40−42 Furthermore, density profiles show that the height of the adsorbed droplets increases from MIL to TIL, Table 2, (elevated by 15 Å from MIL to TIL). Therefore, it can be concluded that going from mono- to multicationic ILs leads to a decrease in the number of ions interacting with the BN surface and results in adsorbed droplets with higher heights. 3.3. ITIM Analysis. The ITIM as a useful tool for analyzing the surface of disordered phases is able to distinguish the molecules that are right at the interface by providing the list of cations and anions that belong to the interfacial region.43 Table 3 summarizes the surface fraction of cations and anions Xion present in the interfacial region calculated as an average over the last 5 ns of the trajectory. It is obvious that the interfacial

Table 2. Characteristics of the Studied Ionic Liquid Nanodroplets on the Monolayer BN Sheet IL/BN

θ (°)

height (Å)

r (Å)

MIL/BN DIL/BN TIL/BN

73 76 86

45 52 67

38 35 30

area decreases, which could be determined from the r values reported in Table 2. This suggests that by increasing the size of the cation by adding the imidazolium rings, their interaction with the BN surface becomes less favorable. 3.2. Density Profiles. In order to investigate the characteristics of the ionic liquid nanodroplets in the vicinity of the monolayer BN sheet, number density profiles of the cation, anion, and IL along the z coordinate, perpendicular to the sheet were analyzed. According to Figure 5, in all cases, the first sharp peak near the BN surface shows a rather rigid or solid-like structure formation due to the presence of the BN surface. However, the first peak (Figure 5c) for the TIL nanodroplet is notably lower than those of MIL and DIL 13555


Article

The Journal of Physical Chemistry C Table 3. Fraction of Cations and Anions Belonging to the Interfacial and Subinterfacial Layers for the Studied ILs IL

Xcation

Xanion

MIL DIL TIL

53 38 32

47 62 68

layer is more populated with cations in all studied nanodroplets and the cation/anion ratio decreases going from MIL to TIL, which is expectable since generally there are larger numbers of anions than cations in multicationic ILs. The dynamics of the ion exchange between the consecutive layers of the liquid phase has been studied using the survival probability L(t) evaluated for ions in the interfacial layer. L(t) is the probability that an ion that is in the given layer at time t0 remains in this layer uninterruptedly up to time t0 + t. Since the exchange of the molecules between the surface layer and the bulk liquid phase is a process of first-order kinetics, the L0(t) and L(t) functions are of exponential decay. L(t) values for the cations and anions of the interfacial layers of the studied ILs are shown in Figure 6. There is a sharp fall of L(t) functions for a very short period of time followed by a steady fall at larger observation time. The L(t) functions are fitted by the sum of two stretched exponential decays for 5 ns of elapsed time functions as follows44 L(t ) = a1exp( −( t τ1 )β1 ) + a 2 exp( −( t τ2 )β2 )

(1)

where τ1 and τ2 are the residence times associated with two exponential functions, a1 and a2 are the contributions from the individual terms, and the β terms stand for stretched exponents. The average residence time ⟨τs⟩ can be expressed as the weighted average of τ1 and τ2 as follows: ⟨τ s⟩ = a1τ1 + a 2τ2

(2)

The parameters of fitting are listed in Table 4. The residence time gives an idea about the mean lifetime of an ion in the interfacial layer, and the difference in the residence time of the ions is related to their different affinities to the surface.45 According to the Table 4, the average lifetime of cations in the interfacial layer is slightly larger than those of anions in all studied nanodroplets, which indicates the higher tendency of cations to occupy the BN surface. Furthermore, the average lifetime of the cation in the interfacial layer decreases from MIL to TIL. The average residence time for the monocation is approximately two times larger than that for the trication, which could be related to the higher affinity of the monocation to the BN surface compared to that of the trication. Furthermore, the very close values of residence times for the anion and cation in the TIL nanodroplet can be attributed to the decreased competition between the cation and anion for occupying the BN surface. 3.4. Orientation of Cations. Density profiles showed that the adsorbed layer on the BN sheet results in a remarkable structuring of the ions in this region, which is different to that in the inner regions of the droplets. In the inner regions of the droplets, the arrangements of the ions are mainly controlled by Coulombic interactions, while in the interfacial region, a competition between the interionic and ionic liquid−nanosheet interactions is responsible for the arrangement of ions at the interface. This arrangement has been quantified through the bivariate orientational method.46 According to the method,

Figure 6. Survival probability of (a) MIL, (b) DIL, and (c) TIL in the interfacial layer for anions and cations of the studied ILs.

two independent spherical polar angles (ϕ and θ) are defined as shown in Figure 7. θ is the angle between the z axis normal to the interface and the bisector vector of the imidazolium ring, rbisector. ϕ is the angle between the imidazolium ring plane normal vector and the projection of z onto the plane made by the two vectors, the normal vector of the imidazolium ring plane surface (n1), and the normal vector of the rbisector−n1 plane (n2). The maps give a contour of ϕ versus cos(θ). 13556


Article


just a limited number of imidazolium rings (mostly middle rings) of the interfacial cations lie parallel to the surface and a large number of the rings (side and middle rings) tend to tumble sideways. This arrangement of ions at the interface for TIL could be explained considering that TIL has bulky cations with large steric effects due to the existence of three cationic moieties in a cation. This complex structure does not allow all three imidazolium rings of cations to adopt a parallel orientation with respect to the BN surface since this orientation can strongly reduce the anion−cation association in the trication. It has been found in our recent study that trications tend to form a helical structure in their bulk phase in order to maximize the hydrogen and electrostatic interactions with anions and minimize the repulsive interactions arising from neighboring cationic moieties.47 Existence of helical structures of the trication in the vicinity of the BN surface is also possible; however, we cannot propose a specific conformation for the trication because there are many other conformational possibilities for such a complex cation in the vicinity of the BN surface. Overall, the results obtained from the ITIM and orientational analysis indicate that while there is a relatively strong adsorption of monocations on the BN monolayer, which is facilitated by the ion packing (monocations can relatively easily tilt their tail for anion packing), the large cation size and steric effects in multications, especially the trication, restricts the parallel ordering of the rings and therefore reduce the packed ion number and cause weaker adsorption on the BN surface. It should be noted that DIL nanodroplets exhibit a denser adsorbed region and stronger adsorption on the BN surface compared to TIL nanodroplets, which is due to the special bending of dications in the vicinity of the BN surface. 3.5. Work of Adhesion and Binding Energy. The percentage of the binding energy of the studied IL nanodroplet to the BN surface, %Ebind, is defined as

Table 4. Parameters of the Stretched Exponential Fit to the Average Survival Time Correlation Function for Cations and Anions of the Studied ILs IL

a1

τ1 (ns)

β1

a2

τ2 (ns)

β2

⟨τs⟩ (ns)

cation (DIL) anion (DIL) cation (MIL) anion (MIL) cation (TIL) anion (TIL)

0.41 0.67 0.47 0.55 0.31 0.28

0.06 0.08 0.09 0.08 0.08 0.07

0.59 0.65 0.77 0.77 0.94 0.92

0.56 0.30 0.40 0.27 0.60 0.62

0.67 0.71 0.65 0.77 0.34 0.35

1.07 1.40 1.29 1.58 1.69 1.54

0.31 0.27 0.30 0.25 0.15 0.14

Figure 7. Molecular fixed coordinates (rbisector, n1, and n2) defined for bivariate orientational analysis and spherical polar angles (ϕ and θ) applied to the ring plane of MIL, DIL, and TIL.

Figure 8a illustrates that for MIL, the most probable region is centered at 0.02 < cos θ < 0.12 with 80 < ϕ < 90°. In this orientation, the plane of the imidazolium ring of IL lies parallel to the BN sheet with a tilted angle of 10°. Parallel stacking of the imidazolium rings at the BN surface is probably due to the π−π stacking interactions between the imidazolium ring and the BN surface. This observation resembles the results computationally obtained for imidazole cations and the benzoate anion in contact with both graphene and BN sheets.31,40 Interestingly, the orientational ordering in DIL is quite different from that of MIL. As shown in Figure 8b, two different orientational preferences can be identified for the imidazolium ring plane of the cation. The most probable regions are centered at −0.1 < cos θ < 0.1 with 40 < ϕ < 60° and −0.1 < cos θ < 0.2 with 75 < ϕ < 90°. Thus, while the bisector in DIL takes the maximum angle of 95° with respect to the z axis, the ring plane tumbles sideways and adopts various orientations. Therefore, a probable orientation for DIL at the interface is that an imidazolium ring lies almost parallel to the BN surface while another ring takes a tilt angle with respect to the surface to leave a space for the anion to maximize both interionic and IL−BN sheet interactions. Previous MD simulations of bulk properties of ILs have shown that the anions form a cluster around the rings of dications and trications.24,25 For TIL, we dissected the cation to middle rings and side rings for orientational analysis because of the structural complexity of the trication. For the middle ring, the most probable regions are centered at −0.3 < cos θ < −0.1 with 60 < ϕ < 75° and 0 < cos θ < 0.2 with 70 < ϕ < 90°, Figure 8c. It can be clearly seen that the orientational ordering in DIL and TIL is lower than that of MIL. Also, for the imidazolium ring plane of the side rings, two different orientational preferences can be identified. The most probable orientation is at −0.6 < cos θ < −0.5 with 35 < ϕ < 50° and 0.8 < cos θ < 0.9 with 0 < ϕ < 10°, Figure 8d. Accordingly, this orientation indicates that

%E bind = 100 ×

[ECFG(BN + drop) − (ECFG(BN) + ECFG(drop)] [(ECFG(BN) + ECFG(drop)] (3)

where ECFG refers to the configurational energy (total potential energy). In all MD simulations, the surface atoms were taken as frozen, so ECFG of the BN sheet does not contribute to the binding energy. Figure 9 compares the percentage of binding energy of the studied nanodroplets on the BN surface calculated from the ensemble average of the last 2 ns of the simulation time. It shows that the binding of MIL nanodroplets with BN surfaces is stronger than those of DIL and TIL nanodroplets. Furthermore, the percentage binding energy decreases by increasing the number of positively charged moieties on the cation from MIL to TIL. Nevertheless, the behavior of %Ebind with regard to the studied ILs is nonlinear. This result is in good consistency with the trend of calculated contact angles for the studied ILs. The work of adhesion (Wadh), which is defined as the reversible free energy change for making free surfaces from interfaces, was obtained by comparing the energy of the relaxed interface with the energy of the relaxed droplets and hBN at infinite separation. Therefore, the work of adhesion can be calculated as follows48 13557


Article


Figure 8. Bivariate orientation distribution of imidazolium rings in the cationic moiety of (a) MIL, (b) DIL, (c) TIL (central ring), and (d) TIL (side rings). Red corresponds to a high normalized probability, and blue corresponds to a low probability.

simulations have also revealed that stronger IL−surface interactions result in the greater wetting of IL nanodroplets and the smaller contact angles.49 It can be concluded that the charge and size of the cationic moiety have an important role in determining the strength of the interaction between the IL and BN surface.

Wadh = [ECFG(BN + drop) − (ECFG(BN) + ECFG(drop)] /A

(4)

where A is the area of the nanodroplets in contact with the BN surface. The calculated work of adhesions for MIL, DIL, and TIL nanodroplets on the BN surface is 101.01, 95.15, and 75.09 mJ/m2, respectively, which shows that by increasing the number of cationic groups from MIL to TIL, the affinity of IL to the BN surface decreases. Therefore, different affinities of the studied ionic liquids to the BN surface could be the origin of the difference in the contact angles. Previous MD

4. CONCLUSIONS Classical MD simulations have been carried out on ionic liquid nanodroplets of three different categories of ILs, that is, MIL, DIL and TIL, in contact with the BN surface. The effect of 13558


Article

The Journal of Physical Chemistry C ORCID

Elaheh Sedghamiz: 0000-0003-3998-9701 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors are indebted to the Institute for Studies in Theoretical Physics and Mathematics (IPM) for providing us the HPC Cluster to perform part of simulations.

■

(1) Liu, Q.; Baoxing, X. Actuating Water Droplets on Graphene via Surface Wettability Gradients. Langmuir 2015, 31, 9070−9075. (2) Wen, L.; Tian, Y.; Jiang, L. Bioinspired Super-Wettability from Fundamental Research to Practical Applications. Angew. Chem., Int. Ed. 2015, 54, 3387−3399. (3) Castejón, H. J.; Wynn, T. J.; Marcin, Z. M. Wetting and Tribological Properties of Ionic Liquids. J. Phys. Chem. B 2014, 118, 3661−3668. (4) Villa, F.; Marengo, M.; De Coninck, J. A new model to predict the influence of surface temperature on contact angle. Sci. Rep. 2018, 8, 6549. (5) Burt, R.; Birkett, G.; Salanne, M.; Zhao, X. S. Molecular dynamics simulations of the influence of drop size and surface potential on the contact angle of ionic-liquid droplets. J. Phys. Chem. C 2016, 120, 15244−50. (6) Š ikalo, Š .; Wilhelm, H. D.; Roisman, I.V.; Jakirlić, S.; Tropea, C. Dynamic contact angle of spreading droplets: Experiments and simulations. Phys. Fluids 2005, 17, No. 062103. (7) Taherian, F.; Marcon, V.; van der Vegt, N. F. A.; Leroy, F. What Is the Contact Angle of Water on Graphene. Langmuir 2013, 29, 1457−1465. (8) Jensen, T. R.; Jensen, M. Ø.; Reitzel, N.; Peters, G. H.; Kjaer, K.; Bjørnholm, T. Water in Contact with Extended Hydrophobic Surfaces: Direct Evidence of Weak Dewetting. Phys. Rev. Lett. 2003, 90, No. 086101. (9) Miwa, M.; Nakajima, A.; Fujishima, A.; Hashimoto, K.; Watanabe, T. Effects of the Surface Roughness on Sliding Angles of Water Droplets on Superhydrophobic Surfaces. Langmuir 2000, 16, 5754−5760. (10) Adamson, A. W.; Gast, A. P. Physical Chemisry of Surfaces, 6th ed.; John Wiley & Sons: New York, 1997 (11) Rahman, A.; Stillinger, F. H. Molecular Dynamics Study of Liquid Water. J. Chem. Phys. 1971, 55, 3336−3359. (12) Li, H.; Zeng, X. C. Wetting and Interfacial Properties of Water Nanodroplets in Contact with Graphene and Monolayer Boron− Nitride Sheets. ACS Nano 2012, 6, 2401−2409. (13) Wei, N.; Lv, C.; Xu, Z. Wetting of Graphene Oxide: A Molecular Dynamics Study. Langmuir 2014, 30, 3572−3578. (14) Weingärtner, H. Angew. Chem., Int. Ed. 2008, 47, 654. (15) Malali, S.; Foroutan, M. Study of Wetting Behavior of BMIM +/PF6−Ionic Liquid on TiO2 (110) Surface by Molecular Dynamics Simulation. J. Phys. Chem. C 2017, 121, 11226−33. (16) Herrera, C.; García, G.; Atilhan, M.; Aparicio, S. Nanowetting of graphene by ionic liquid droplets. J. Phys. Chem. C 2015, 119, 24529−37. (17) Li, G.-X.; Liu, Y.; Wang, B.; Song, X.-M.; Li, E.; Yan, H. Preparation of Transparent BN Films with Superhydrophobic Surface. Appl. Surf. Sci. 2008, 254, 5299−5303. (18) Lee, C. H.; Drelich, J.; Yap, Y. K. Superhydrophobicity of Boron Nitride Nanotubes Grown on Silicon Substrates. Langmuir 2009, 25, 4853−4860. (19) Garcia, G.; Atilhan, M.; Aparicio, S. Adsorption of choline benzoate ionic liquid on graphene, silicene, germanene and boronnitride nanosheets: a DFT perspective. Phys. Chem. Chem. Phys. 2015, 17, 16315−16326.

Figure 9. IL−BN surface interaction energy percentage calculated for MIL, DIL, and TIL nanodroplets on the BN surface. Lines are reported for guiding purposes.

increasing number of positively charged centers of the cation on the wettability and interfacial properties of ILs is analyzed, showing that this increase leads to remarkable changes in wettability and structuring of ILs on the BN surface. Calculated contact angles of all ionic liquid nanodroplets on the BN surface demonstrate favorable wettability of the studied ILs on the BN sheet and an increase going from MIL to TIL. Density profiles show that placement of the MIL nanodroplet on the BN surface occurs with a dense layer of cations and anions. In addition, they demonstrate that although both cations and anions exist at the interface, the cations lie at closer distances to the BN surface. ITIM results show that the average survival time of cations in the vicinity of the surface layer is larger than that of anions. Furthermore, the probability for the cations to be found in the interfacial layer at time t0 and t0 + t shows a decrease in the average survival time with increasing number of imidazolium moieties of the cations. The bivariate orientational analysis reveals that the imidazolium rings in MIL tend to lie flat at the BN surface, while this orientational preference decreases going from MIL to TIL. A weak layering of trications on the BN surface is observed compared to those of the monocation and dication, which is attributed to special steric effects caused by existence of three positively charged centers on the trication. Furthermore, percentage binding energies and adhesion work decreases by increasing the number of positively charged centers of the cation from MIL to TIL. It can be concluded that different affinities of the studied ionic liquids to the BN surface could be the origin of the difference in the contact angles. These observations represent the important role of the number of charged moieties of the cation, steric effects, and structural restrictions when the adsorption and wetting of ILs are concerned.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.8b11987. Tables of electrostatic charges for the atoms of studied MIL, DIL, and TIL and Lennard-Jones parameters and figure of error bars of contact angles (PDF)

■

REFERENCES

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. 13559


Article

The Journal of Physical Chemistry C (20) Shakourian-Fard, M.; Kamath, G.; Jamshidi, Z. Trends in Physisorption of Ionic Liquids on Boron-Nitride Sheets. J. Phys. Chem. C 2014, 118, 26003−26016. (21) Sharma, P. S.; Payagala, T.; Wanigasekara, E.; Wijeratne, A. B.; Huang, J.; Armstrong, D. W. Trigonal Tricationic Ionic Liquids: Molecular Engineering of Trications to Control Physicochemical Properties. Chem. Mater. 2008, 20, 4182−4184. (22) Shirota, H.; Mandai, T.; Fukazawa, H.; Kato, T. Comparison between Dicationic and Monocationic Ionic Liquids: Liquid Density, Thermal Properties, Surface Tension, and Shear Viscosity. J. Chem. Eng. Data 2011, 56, 2453−2459. (23) Anderson, J. L.; Ding, R.; Ellern, A.; Armstrong, D. W. Structure and Properties of High Stability Geminal Dicationic Ionic Liquids. J. Am. Chem. Soc. 2005, 127, 593−604. (24) Sedghamiz, E.; Moosavi, M. Tricationic Ionic Liquids: Structural and Dynamical Properties via Molecular Dynamics Simulations. J. Phys. Chem. B 2017, 121, 1877−1892. (25) Moosavi, M.; Khashei, F.; Sedghamiz, E. Molecular dynamics simulation of geminal dicationic ionic liquids [Cn(mim)2][NTf2]2 − structural and dynamical properties. Phys. Chem. Chem. Phys. 2018, 20, 435. (26) Sedghamiz, E.; Moosavi, M. Probing the tricationic ionic liquid/vacuum interface: insights from molecular dynamics simulations. Phys.Chem.Chem.Phys. 2018, 20, 14251−14263. (27) Smith, W.; Todorov, I. DL_POLY Classic Molecular Simulation Package; Daresbury Laboratory: Cheshire, U.K., 2012. (28) Humphrey, W.; Dalke, A.; Schulten, K. VMD - Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (29) Canongia Lopes, J. N.; Pádua, A. A. Molecular Force Field for Ionic Liquids III: Imidazolium, Pyridinium, and Phosphonium Cations; Chloride, Bromide, and Dicyanamide Anions. J. Phys. Chem. B 2006, 110, 19586−19592. (30) Tsuzuki, S.; Shinoda, W.; Saito, H.; Mikami, M.; Tokuda, H.; Watanabe, M. Molecular Dynamics Simulations of Ionic Liquids: Cation and Anion Dependence of Self-Diffusion Coefficients of Ions. J. Phys. Chem. B 2009, 113, 10641−10649. (31) Torkzadeh, M.; Moosavi, M. A combined molecular dynamics simulation and quantum mechanics study on the physisorption of biodegradable CBNAILs onh-BN nanosheets. The J. Chem. Phys. 2018, 149, No. 074704. (32) Albe, K.; Müller, W. Modelling of boron nitride: Atomic scale simulations on thin film growth. Comput. Mater. Sci. 1998, 10, 111− 115. (33) Kamath, G.; Baker, G. A. Are ionic liquids suitable media for boron nitride exfoliation and dispersion? insight via Molecular Dynamics. RSC Adv. 2013, 3, 8197−8202. (34) Won, C. Y.; Aluru, N. R. Structure and Dynamics of Water Confined in a Boron Nitride Nanotube. J. Phys. Chem. C 2008, 112, 1812−1818. (35) Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511−519. (36) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, U.K., 1987 (37) Darvas, M.; Pojják, K.; Horvai, G.; Jedlovszky, P. Molecular dynamics simulation and identification of the truly interfacial molecules (ITIM) analysis of the liquid-vapor interface of dimethyl sulfoxide. J. Chem. Phys. 2010, 132, 134701. (38) Hantal, G.; Terleczky, P.; Horvai, G.; Nyulászi, L.; Jedlovszky, P. Molecular Level Properties of the Water−Dichloromethane Liquid/Liquid Interface, as Seen from Molecular Dynamics Simulation and Identification of Truly Interfacial Molecules Analysis. J. Phys. Chem. C 2009, 113, 19263−19276. (39) Giovambattista, N.; Debenedetti, P. G.; Rossky, P. J. Effect of Surface Polarity on Water Contact Angle and Interfacial Hydration Structure. J. Phys. Chem. B 2007, 111, 9581−9587. (40) Wang, S.; Li, S.; Cao, Z.; Yan, T. Molecular Dynamic Simulations of Ionic Liquids at Graphite Surface. J. Phys. Chem. C 2009, 114, 990−995.

(41) Sieffert, N.; Wipff, G. Solvation of Sodium Chloride in the 1Butyl-3-methyl-imidazolium Bis(trifluoromethylsulfonyl)imide Ionic Liquid: A Molecular Dynamics Study. J. Phys. Chem. B 2007, 111, 7253−7266. (42) Sha, M.; Zhang, F.; Wu, G.; Fang, H.; Wang, C.; Chen, S.; Zhang, Y.; Hu, J. Ordering layers of [Bmim][PF6] ionic liquid on graphite surfaces: molecular dynamics simulation. J. Chem. Phys. 2008, 128, 134504. (43) Palchowdhury, S.; Bhargava, B. L. Surface Structure and Dynamics of Ions at the Liquid−Vapor Interface of Binary Ionic Liquid Mixtures: Molecular Dynamics Studies. J. Phys. Chem. C 2016, 120, 5430−5441. (44) Pártay, L. B.; Hantal, G.; Jedlovszky, P.; Vincze, Á .; Horvai, G. A new method for determining the interfacial molecules and characterizing the surface roughness in computer simulations. Application to the liquid−vapor interface of water. J. Comput. Chem. 2008, 945−956. (45) Idrissi, A.; Hantalbc, G.; Jedlovszky, P. Properties of the liquid−vapor interface of acetone−methanol mixtures, as seen from computer simulation and ITIM surface analysis. Phys. Chem. Chem. Phys. 2015, 17, 8913−8926. (46) Pártay, L.; Jedlovszky, P.; Jancsó, G. Orientation of the 3methylpyridine molecules at the liquid−vapour interface of their aqueous solution. Chem. Phys. Lett. 2006, 420, 367−371. (47) Sedghamiz, E.; Khashei, F.; Moosavi, M. Linear tricationic ionic liquids: Insights into the structural features using DFT and molecular dynamics simulation. J. Mol. Liq. 2018, 271, 96−104. (48) Nikkhah, S. J.; Moghbeli, M. R.; Hashemianzadeh, S. M. A molecular simulation study on the adhesion behavior of a functionalized polyethylene-functionalized graphene interface. Phys. Chem. Chem. Phys. 2015, 17, 27414−27. (49) Guan, Y.; Shao, Q.; Chen, W.; Liu, S.; Zhang, X.; Deng, Y. Dynamic Three-Dimensional Nanowetting Behavior of ImidazoliumBased Ionic Liquids Probed by Molecular Dynamics Simulation. J. Phys.Chem. C. 2017, 121, 23716−23726.

13560


Molecular Dynamics Simulation of Wetting and ... - ACS Publications

Recommend Documents