Article Cite This: Langmuir XXXX, XXX, XXX-XXX
pubs.acs.org/Langmuir
Toward High-Throughput Computational Screening of Carbon Nanotube Solvents A. Hardy and H. Bock* Institute of Chemical Sciences, Heriot Watt University, Edinburgh EH14 4AS, U.K. S Supporting Information *
ABSTRACT: We use the corresponding distances method (CDM) to computationally assess the quality of 10 experimentally tested carbon nanotube (CNT) solvents. The CDM produces accurate and high-resolution potential of mean force curves from a single simulation per solvent. The method’s very high efficiency allows us to investigate an unprecedented number of solvents in one study. The simulation results indicate that none of the tested molecules are solvents in the thermodynamic sense; instead, they are dispersants preventing reaggregation of already dispersed CNTs. We find that the dispersion free energy barrier correlates very well with the experimentally measured performance of the dispersants; i.e., the simulations place the solvents in the correct performance order. Our analysis of the structure−function relationship rationalizes this order.
■
INTRODUCTION Carbon nanotubes (CNTs) have a range of properties that make them potentially useful for a number of applications: they are one of the strongest known materials (in terms of tensile strength and elastic modulus1), their electronic properties vary between metallic and semiconducting, and they have excellent thermal properties. Despite being known for much more than two decades,2 they have struggled to find their way into real world, commercial applications. The key bottleneck remains the lack of an efficient method for turning the intractable black raw material, consisting of bundled and entangled tubes, into a dispersion or, better, solution of individual tubes. The strong aggregation tendency of CNTs has its origin in the very strong tube−tube cohesion of approximately −40kT nm−1 for a pair of (10,10) tubes (Ø = 1.36 nm).3 As carbon nanotubes can be centimeters long,4 this cohesive force can accumulate to an extremely large attraction, making exfoliation, so far, an insurmountable challenge. Unfortunately, CNTs realize their full potential only when they are dispersed into individualized tubes. The dispersion problem has been tackled using a range of liquid phase processing techniques,5−7 but despite years of work, high-yield individualization is still the most significant barrier to the widespread adoption of CNTs as a commercial material. A surfactant in water is the most common method used to create CNT dispersions. Water is environmentally benign, and surfactants have historically been successful in creating dispersions of spherical colloids.8 Some popular surfactants for use with CNTs include sodium dodecyl sulfate (SDS),9,10 sodium dodecyl benzenesulfonate (SDBS),9−11 and Triton X100.10,12 Despite being in use for more than a decade, these © XXXX American Chemical Society
surfactant-based dispersion systems still suffer from poor yields, often well under 10% of the initial CNT mass, with CNT concentrations in the 0(1) mg/mL range.9−11,13 Computer simulations have also shown that surfactants that form spherical aggregates on the tubes, as SDS does, are ineffective dispersants.14,15 Surfactants also introduce processing problems: they leave behind a residue that is difficult to remove and can be detrimental to the properties of the final materials.16 Pure solvents avoid this problem, as long as the molecules are not too large. Solvents have been studied extensively.17−21 We focus particularly on work by Bergin et al., 18 who experimentally tested dozens of solvents, repeatedly measuring the amount of dispersed tubes for each. One of the key results of this work was finding one particular molecule, CHP (cyclohexyl-pyrrolidinone), outperforming all other solvents by a factor of 5, dispersing 3.5 mg of CNTs in 1 mL of solvent. Understanding why CHP performs so much better than all other solvents would be an important step toward the rational design of solvents that achieve even higher CNT loadings. Fundamentally, there are two ways to stabilize individualized carbon nanotubes in a liquid: dispersion and spontaneous dissolution. In a true thermodynamic solvent, the free energy of the individualized tubes immersed in the solvent is lower than that of the bundle. Thus, the solution is thermodynamically stable. In a dispersant, the free energy of the aggregated state (CNT bundle) is lower than that of the dispersion. The tubes Received: July 25, 2017 Revised: September 27, 2017 Published: October 2, 2017 A
DOI: 10.1021/acs.langmuir.7b02600 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir can, however, be kept in the metastable dispersion if the dispersant creates a sufficiently high free energy barrier between the dispersed and aggregated states. The reaggregation barrier needs to be large enough that it cannot be overcome by the system’s thermal energy. The reaggregation barrier necessarily also creates a barrier against dispersing the tubes from the bundle. This is inconvenient, but the proposed unzippering mechanism22 would provide a route from the bundled state to the individualized state that would essentially bypass the dispersion barrier. Only the initialization of the process requires the crossing of a barrier. Current experimental methods almost exclusively use ultrasonication for this, but the applied energy densities are so high that it breaks carbon−carbon bonds within the tubes, thereby degrading their unique properties.21,23 Unfortunately, the principle of the unzippering mechanism also applies to the reaggregation process. This has the important consequence that if a tube has not been completely exfoliated from a bundle, the exfoliated part cannot be kept in the metastable state by any dispersant and will always spontaneously rejoin the bundle. These problems motivate the search for a true thermodynamic solvent, in which the individualized tubes are thermodynamically stable. In a thermodynamic solvent, the exfoliation process would proceeded spontaneously via the unzippering mechanism or lateral sliding once the respective initial barrier has been overcome. In computer simulations, free energy differences and barrier heights are readily available via the potential of mean force. Here we use the corresponding distances method (CDM) to asses the solvent and dispersant quality of a number of solvents that have previously been tested experimentally.18 In particular, we investigate if there is any indication that these solvents could be solvents in the thermodynamic sense. We include a number of (experimentally) poor and the very best solvents to elucidate the molecular properties that are responsible for the difference.
■
Figure 1. Schematic of a typical simulation system setup (top) and the CDM system (bottom). The top panel is provided only to demonstrate the correspondence of the CDM to the parallel tube simulations and does not represent any system we actually simulated. The two systems are linked via the corresponding distances d = l + d′, where d′ is the distance between the tube surface and the first solvent layer.
Figure 2. 10,10 nanotube highlighting one unit cell. There are two rings of C atoms that are differentiated only by their relative rotation around the cylinder axis.
MODEL AND SIMULATION
We use the CDM24,25 to calculate the proton motive force (PMF) of parallel 10,10 carbon nanotubes for a range of solvents, as it produces highly accurate and high-resolution PMF curves from only a single simulation per solvent. The PMF w(d) is obtained from mean force F(d) on the nanotubes at tube−tube distance d via integration
w(d) = −
∫d
a diameter of 1.36 nm. The nanotubes were built with the aid of the Buildcstruct script.37 The nanotubes are approximated as rigid rods. To achieve this, the CNT atoms are “frozen”; i.e., forces on the atoms are calculated as normal, but their coordinates are not updated. The CNT atoms were given the nonbonded characteristics of a “bare carbon” atom in the GROMOS 53a6 force field. GROMOS 53a6 is shown to reproduce the tube−tube potential for (10,10) carbon nanotubes well.3 The initial simulation box dimensions were 10.739 nm × 5 nm × 122.63 nm, containing two 10,10 single-walled carbon nanotubes. The axis of the first tube goes through the box center and is oriented parallel to the z-axis, whereas the second is rotated around y-axis from the same position by 5°. This small angle guarantees the very high accuracy of the CDM approximation and provides a very high data density of F(d) of 400 data points per σ of tube−tube distance, which is vital for accurate numerical integration of the mean force to obtain the PMF. The longer, diagonal tube was 123.10 nm long, and the slightly shorter horizontal tube was 122.97 nm long. At the desired box size for our simulations, it was not possible to simultaneously maintain the periodic structure of both tubes under periodic boundary conditions, and small artifacts arose near the boundaries. This has been accounted for by removing 50 rows of carbon atoms either side of the boundary. The simulation box is filled with solvent by overlaying pre-equilibrated solvent boxes and removing any overlapping solvent molecules. Any solvent molecule inside the nanotubes is also removed, so the tubes are simulated without any solvent inside.
∞
F(d) ∂d
(1)
In this work, d is always the surface-to-surface distance (top panel of Figure 1) and F(d) and d have the same direction; i.e., a positive F points in the direction of an increasing d. In practice, the forces are integrated from a “far enough” distance, d*, chosen so that F(d) = 0 beyond d*.
w(d) = −
∫d
d*
F(d) ∂d
(2)
The CDM for atomistic simulations has been described in detail in ref 24. In brief, instead of aligning the tubes in parallel, they form an acute angle and cross through each other in the center of the simulation box (Figure 1). Then by extraction of the local forces for each ring of carbon atoms (Figure 2), the entire mean force curve is obtained. Molecular dynamics simulations were performed using the GROMACS 5.1 simulation package.26−32 The force field used was GROMOS 53a6.33 Parametrization of the solvents was performed with the aid of the Automated Topology Builder34−36 and manually verified. The tubes used in this study are (10,10) armchair CNTs with B
DOI: 10.1021/acs.langmuir.7b02600 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir The system is initially energy minimized, followed by a temperature equilibration phase in the NVT ensemble and an equilibration phase in the NPysxszT ensemble, where the box dimensions in the plane of the tubes (x and z) are fixed, and the out-of-plane dimension (y) is coupled to the y-component of the pressure. The average box dimensions of the last 1 ns of the pressure-coupled equilibration were used for the final data collection phase in the NVT ensemble (of which the first 1 ns of data is discarded to allow for equilibration). The simulation parameters for each step can be found in Table 1.
Table 1. Simulation Parametersa production ensemble simulation time (ns) time step (ps) integrator van der Waals and electrostatic cutoff (nm) van der Waals scheme electrostatic scheme thermostat thermostat coupling time/ constant (ps) van de Waals mixing rule machine precision neighbor list update frequency (steps) neighbor list cutoff (nm) isothermal/isobaric equilibration ensemble (see text) simulation time (ns) barostat barostat coupling constant (ps) isothermal compressibility (bar−1) isothermal equilibration time (ns) ensemble energy minimization EM type maximum number of steps Fmax (kJ mol−1)
NVT 10 0.002 velocity verlet 1.6
Figure 3. Ten solvents tested. Further information on each can be found in Table S1.
Lennard-Jones 12−6, potential shift reaction field zero, potential shift verlet velocity rescale 0.1
We begin by discussing a single solvent to introduce the relevant quantities needed for the comparison. For this purpose, we choose CBrCl3, which is discussed in great detail in ref 24. In Figure 4a, three PMF curves are shown: the direct tube−tube PMF in vacuum (ttPMF) in blue, the solventmediated PMF (smPMF) in magenta, and the sum of the two, i.e., the actual PMF, in yellow. The horizontal axis represents
geometric mean double 20 1.9 NPzsxyT 2 Berendsen 2.0 4 × 105 0.3 NVT conjugate gradient 10000 103
a
Parameters for the energy minimization and equilibration steps are the same as for the production step unless otherwise stated. Descriptions of the simulation parameters (integrator, nonbonded schemes, etc.) can be found in the GROMACS manual.38
All four CNT branches are used for data collection. It should be noted that these branches may not be identical depending on where exactly the pivot point is placed. It should also be noted that because the repeat unit of the tubes contains two rings of carbon atoms with different atom positions around the tubes (Figure 2), the raw mean force has high-frequency oscillations, whereas for each ring type separately, it is very smooth. This leads to eight separate force curves that we integrate separately. Finally, we average the four contributions for each ring type and then add the two contributions to obtain the PMF.
■
Figure 4. (a) PMF curves for a pair of 10,10 CNTs in CBrCl3, including the tube−tube part (ttPMF, blue), the solvent-mediated contribution (smPMF, magenta), and the full PMF (yellow). (b) Expanded (from the dashed area in panel a) and annotated PMF illustrating the key properties of the PMF, including the reaggregation barrier (wreagg), the dispersion barrier (wdisp), and the stability of the aggregate (wagg).
RESULTS AND DISCUSSION Basis of the Comparison. We study a range of solvents shown in Figure 3 and selected to include the best and some of the poorest candidates from the large list of solvents that have previously been experimentally tested by Bergin et al.18 C
DOI: 10.1021/acs.langmuir.7b02600 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir
Figure 5. Atom density maps of CBrCl3 around a pair of tubes showing the layering of chlorine (a, green), bromine (b, red), and carbon (c, blue). While many layers are visible, the first layer is much more dense and stratified, indicating relatively strong adhesion and orientational order.
tube−tube surface-to-surface distance d. The ttPMF has a very deep minimum at d = 0.33 nm demarcating the aggregated state. The dispersed state is, naturally, characterized by large values of d (see Figure 4b). The solvent-mediated PMF oscillates at distances d ≳ 0.8 nm, and this structure is inherited by the combined PMF. These oscillations originate from layering of the solvent around the nanotubes (Figure 5). Between the tubes these layers are confined. At zero force, a given number of solvent layers are stable and the smPMF has a minimum. Increasing d from this point stretches the structure and results in an attractive restoring force. Reducing d instead compresses the layers, resulting in a solvent-mediated repulsion. Further compression leads to complete squeeze-out of the central layer, and a new stable tube−tube distance is found. Because of the high degree of curvature of the tubes, this process is rather gradual. Layers beyond the first are weak and require only a few kT nm−1 to be squeezed out. The first layer is much denser than all the others, and its molecules are orientationally ordered, which can be surmised from the stratification of the Cl and Br densities (Figure 5). The high density suggests strong adhesion to the tubes, and as a consequence, it is much more difficult to squeeze out the first layer. This is evident by the strongly repulsive region in the smPMF at d ≲ 0.8 nm. As this is the region of strong tube−tube cohesion, the strong adherence of the first solvent layer to the tubes is the cause of almost all of the solvent-mediated repulsion that counteracts tube aggregation. Thus, the structure of the first layer is crucial for a solvent’s performance in CNT dissolution and dispersion. Combining the smPMF, with its broad repulsive region at small d, with the deep attractive well of the ttPMF leads to a PMF for tubes immersed in CBrCl3 that has three important features, illustrated in Figure 4b. First, the stability of the aggregated state has been significantly weakened (wagg = 10.4kT nm−1 at d = 0.33 nm), but it remains the thermodynamically stable state, by a large margin. Thus, CBrCl3 is clearly not a thermodynamic solvent. For this, wagg would have to be very close to zero, or, better yet, positive. Second, because the range of the smPMF repulsion is longer than that of the ttPMF attraction, a barrier of wreagg = 12.4kT nm−1 forms in the PMF at d = 0.58 nm that separates the aggregated and dispersed states. Such a reaggregation barrier is the defining feature of a dispersant. It is immediately clear that for parallel approach of two tubes CBrCl3 would form a quite enormous barrier, even if the tubes are only 10 nm long. However, at an angle, the overlap area is small, and therefore, the barrier is equally small and may be overcome by the thermal energy of the system. Thus, CBrCl3 is neither a thermodynamic solvent nor a dispersant for CNTs, which is consistent with the experimental results of Bergin et al.18
The third feature of the PMF that is relevant for the performance of a thermodynamic solvent or dispersant is the dispersion barrier. The dispersion barrier is created by the same maximum in the PMF that leads to the reaggregation barrier; it is only climbed from the other side. For CBrCl3, the dispersion barrier is very large (wdisp = 24.6kT nm−1). This barrier is very high and the reason why ultrasonication is needed to activate the dispersion process. After that, dispersion may proceed spontaneously following the unzipping mechanism,22 but that would require the dispersed state to be thermodynamically stable; otherwise, the reverse, i.e., aggregation, is spontaneous. Thus, good dispersants have the intrinsic problem that they make dispersion more difficult. This highlights again the need to find a thermodynamic solvent for CNTs. Comparison of Solvents and Correlation to Experimental Results. We used experimental data obtained from ref 18 as a comparison to our simulated data. In these experimental solubility tests, a number of solvents, including CHP and DMPU, performed much better than CBrCl3 did. They are therefore included in this study. The PMF for each of the selected solvents is shown in Figure 6. Given the large variation
Figure 6. PMFs of all simulated solvents. Differences are small but large enough to have a strong effect on solvent/dispersant quality (see the text).
in solubility test performance among these solvents, it is surprising that all PMFs are very similar. One may, at first, expect that at least CHP, which achieves a CNT dispersion loading (Cmax) 5-fold higher than that of the next best solvent, should stand out, but it does not. However, as we discuss below, small differences in the PMF lead to large differences in performance. Perhaps the most important observation from Figure 6 is that while all solvents significantly reduce the stability of the aggregated state (from wagg = −44kT nm−1 in vacuum), it D
DOI: 10.1021/acs.langmuir.7b02600 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir always remains much more stable than the dispersion; i.e., none of these “solvents” are solvents in the thermodynamic sense. If these are not thermodynamic solvents, then it follows that they can at best be dispersants. Indeed, all of them generate a reaggregation barrier. However, as a consequence, they also have a dispersion barrier. This begins to explain why exfoliation of CNTs is so difficult: for all tested molecules, the aggregated state is very stable and all of them generate a high dispersion barrier that cannot be thermally activated. Thus, for all of these molecules, activation of the dispersion process requires energy input and exfoliation requires continuous energy input. As no methodology exists to provide this latter energy in the required form, the search for better dispersants and solvents is critical. Upon closer inspection, we see that CHP has the highest reaggregation barrier, consistent with the experimental results. In colloidal suspensions, reaggregation is thermally activated. If we assume an Arrhenius type behavior for this process, we can link lifetime τ of the dispersion to the reaggregation barrier (wreagg) via τ∝
1 ∝ e wreagg / kT kr
(3)
where kr is the rate constant. If we assume that the preexponential factor is not too dissimilar for the different solvents, we can derive an expression for the lifetime ratio for two solvents (A and B) at the same CNT loading. A B τA ∝ e(wreagg − wreagg)/ kT τB ∼
(4)
Thus, the ratio of the dispersion lifetime for two solvents is related to the difference in barrier height via an exponential function that amplifies small differences in barrier height to large differences in dispersion lifetimes. For example, CHP has a reaggregation barrier of 16.6kT nm−1 and CBrCl3 one of 12.4kT nm−1. Crucially, these figures are per nanometer of tube, meaning that the lifetime ratio increases very quickly with contact length. Additionally, the barriers for perpendicular contact, assumed here to be equivalent to that of 1 nm of tube, are right at the limit where a dispersion becomes stable. This suggests much more stable CNT suspensions in CHP than in CBrCl3, consistent with the experimental results. This analysis suggests a correlation of the computed dispersion barriers for all molecules with the respective experimentally measured CNT loadings using a logarithmic scale for the experimental results. From Figure 7a, it is immediately clear that there is a strong correlation between the two sets of data. In particular, the best and therefore most important of the dispersants, at the top right of the plot, are largely sorted into the correct order. It is interesting to check for any correlation of the other two relevant properties of the PMF with the experimental results. The data in Figure 7c demonstrate that for all solvents the aggregated state is very stable and that differences between the solvents are small. There is also no clear correlation with the experimental solubility test results. However, some of the best dispersants also have the least stable aggregates. This is consistent with our understanding that the repulsive shortrange part of the smPMF is responsible for both effects. If the dispersion process is critical for the performance of the solvents, one would expect an inverse correlation between the size of the dispersion barrier and Cmax; i.e., solvents with a
Figure 7. (a) Reaggregation barrier, (b) dispersion barrier, and (c) stability of the aggregated state in all solvents against experimentally achieved CNT loading Cmax on a log scale, clearly showing the excellent correlation between the reaggregation barrier and the experimental results. Solvents with Cmax = 0 are indicated with a light purple square and have been placed (artificially) on the vertical axes. The spreads of data over 10 simulations for CHP and CBrCl3 are colored green and blue, respectively, with the original data point for each colored purple.
smaller dispersion barrier would disperse more tubes. However, in all cases, the barrier is very high, and we, again, do not find a correlation. It may be noted that the solvent with the highest Cmax also has the highest dispersion barrier; this finding contradicts the hypothesis but is consistent with the fact that the dispersion barrier is intimately linked to the other two properties. These results suggest once more that the experimental systems are indeed suspensions rather than solutions. They also suggest that the dispersion barrier computed using the CDM is the correct objective function for dispersant design. Data Quality. Because relatively small differences in the reaggregation barrier height can have a large impact on a dispersant’s predicted performance, it is important to asses the uncertainty of the data. To do this, 10 repeat runs were E
DOI: 10.1021/acs.langmuir.7b02600 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir performed for CHP and CBrCl3, the best and one of the poorest dispersants in the solubility test, respectively. In addition, CHP is one of the largest molecules in our study, whereas CBrCl3 is one of the smallest. Rather than calculating a variance, we show all data points alongside the rest of the study data in Figure 7. For CBrCl3, the reaggregation barrier data spread over ∼0.5kT nm−1, whereas the range is approximately 2kT nm−1 for CHP. The most important conclusion from this is that the uncertainty in the data does not affect the general correlation between the reaggregation barrier height and the experimental solubility test results. The difference in the spread of the data is most likely related to the size of the molecules. Because of its larger size, CHP has orientational and translational relaxation times much longer than those of the smaller CBrCl3, meaning CBrCl3 can sample more of its phase space than CHP within the time frame of a simulation. The spread of the data for the aggregate PMF is also ∼2kT nm−1. This is ∼15% of the absolute value of the aggregate PMF and therefore good enough to make the distinction between a good and a bad thermodynamic solvent. Structure−Function Relationship. We have so far been able to phenomenologically link the amount of dispersed tubes in the experiments to our model systems via the height of the reaggregation barrier in the PMF. The shape of the relevant region of the PMF is determined by the stability of the first layer, and this is in turn determined by three things: the interactions between the molecules in the first layer and the tube surface, the interactions between the orientationally ordered molecules within the first layer, and the interactions of the first layer with the bulk. All of these are ultimately dictated by the chemical structure of the solvent. To advance beyond a trial and error approach to CNT solvent and dispersant design, we need to understand the link between the chemical structure of a solvent molecule and the resulting PMF. Because of the complexity of the structure−function relationship, we have limited our analysis to a subset of solvents, including CBrCl3, DMPU, and CHP (Figure 7b,c,e). The first performed poorly in both experiment and simulation, and the latter two performed well. As the model tubes have no partial charges and are not polarizable, they interact with the solvent only via van der Waals forces, modeled here by a 12−6 Lennard-Jones (LJ) potential. The strength and range of the LJ potential depend on the interacting pair of atoms. A few relevant examples are shown in Figure 8. It should be noted that although the depth of the attractive potential well varies by a factor of >2, its position shifts by
