Tad
Hogg
*^{a},
Matthew S.
Moses
*^{b} and
Damian G.
Allis
*^{c}
^{a}Institute for Molecular Manufacturing, Palo Alto, CA, USA. E-mail: tad@imm.org
^{b}Independent Consultant, Lafayette, CO, USA. E-mail: matt.moses@jhu.edu
^{c}Department of Chemistry, Syracuse University, Syracuse, NY, USA. E-mail: dgallis@syr.edu
First published on 16th May 2017
A computationally-efficient method for evaluating friction in molecular rotary bearings is presented. This method estimates drag from fluctuations in molecular dynamics simulations via the fluctuation–dissipation theorem. This is effective even for simulation times short compared to a bearing's energy damping time and for rotation speeds comparable to or below typical thermal values. We apply this method to two molecular rotary bearings of similar size at 300 K: previously studied nested (9,9)/(14,14) double-walled carbon nanotubes and a hypothetical rotary joint consisting of single acetylenic bonds in a rigid diamondoid housing. The acetylenic joint has a rotational frictional drag coefficient of 2 × 10^{−35} kg m^{2} s^{−1}. The friction for the nested nanotubes is 120 times larger, comparable to values reported by previous studies. This fluctuation-based method could evaluate dissipation in a variety of molecular systems with similarly rigid and symmetric bearings.
Design, System, ApplicationReducing friction is important when designing molecular machines. One approach uses computer simulations to estimate friction in various machines and identify the best candidates to synthesize. These simulations start a proposed machine in motion and measure how rapidly it slows. Unfortunately, accurate measurements can require many long simulations because of thermal fluctuations in the speed. However, friction ultimately arises from those fluctuations. So instead of measuring friction directly, a more efficient method infers the friction from the fluctuations. We illustrate this method for two molecular rotary bearings: nested nanotubes and a new bearing based on rotation around single carbon–carbon bonds. With no sliding surfaces, the single-bond bearing has about 100 times less friction than nested nanotubes. Thus this bearing could greatly reduce energy dissipation in future molecular machines. |
A key performance measure for molecular machines is their energy dissipation. Unfortunately, experimental measurements are often difficult to obtain, and impossible for molecular machines that cannot yet be fabricated. Computational studies are a reasonable alternative for evaluating and optimizing designs for various applications. Conceptually, estimating dissipation for a molecular machine is straightforward: simulate the machine's behavior over enough time for energy introduced into the machine's operational degrees of freedom to irreversibly transfer to many other degrees of freedom. Unfortunately, such computational studies are challenging for machines with many atoms and when the machine's operation is only weakly coupled to other atoms, leading to low friction and long damping times. Simulating such machines long enough to observe this damping may not be feasible, even for machines operating in vacuum, where the absence of gas or solvent molecules greatly reduces the computational cost. Moreover, even if simulating a single design is feasible, identifying low-power designs may require examining many machine variants, leading to a computationally costly overall design task.
One approach to reduce the computation time is to simulate the machine at speeds well above the typical thermal speed of the machine's motion. In this case, thermal fluctuations are relatively small, allowing definite estimates of damping directly from the change in speed, even over times short compared to the damping time. For example, this approach has been applied to determine friction in nested nanotubes^{21} However, this method may not apply to large molecular machines intended to operate near or below thermal speeds. This is because the machine may behave very differently at low and high speeds. For instance, machines at high speeds could involve forces strong enough to significantly distort components from their low-speed geometries or even break chemical bonds. More generally, high-speed operation may excite high-frequency vibrational modes that would not be relevant at intended operational speeds, thereby significantly overestimating the damping the machine would experience at its intended operating speed.
To address this challenge, this paper examines the effectiveness of an alternate method:^{22} using fluctuations in the simulated behavior to infer friction via the fluctuation–dissipation theorem. In contrast to measuring the damping directly from the decreasing energy during the simulation, fluctuations are evident even in simulations over times much shorter than the machine's damping time and when operating at speeds comparable to or less than thermal speeds.
This paper evaluates this fluctuation-based method for two examples. The first, rotating nested nanotubes, has been studied extensively and allows comparison with direct estimation of dissipation from the decrease of the energy. In this case, the relatively small number of atoms and high dissipation make both the direct and fluctuation-based methods feasible with modest computational cost.
The second example is a hypothetical molecular rotor using acetylenic bonds. These bonds give exceptionally low rotary friction, and thereby would require simulating the machines over long times to reliably measure dissipation directly. Moreover, including enough atoms around those bonds to capture interactions with the rotor's housing results in a significantly larger number of atoms than used for the nanotube example. Thus this rotor is an example where the fluctuation-based method could provide a significant benefit.
Specifically, after discussing related studies, section 3 describes the fluctuation-based method. We then present the two examples. The concluding section discusses future directions for applying the method and its possible limitations. The appendix provides details of the simulations.
In ref. 21, four different MD simulation configurations using the LAMMPS code and AIREBO force field were used to obtain consistent friction estimates for a rotating nested nanotube bearing. Other dissipation estimates for rotating nested nanotube bearings have been determined using intralayer interactions based on the Brenner potential^{23} with interlayer interactions based on the Kolmogorov–Crespi registry-dependent potential,^{24} and with the COMPASS force field.^{25} Dissipation estimates for linear sliding nested nanotube bearings have been determined using a custom force field and custom numerical simulator in ref. 26 and 27, and a custom force field based on ref. 28 and the r-RESPA integrator in ref. 29.
A MD model of experimentally realized graphene-on-graphene sliding is presented in ref. 15. The self-retracting motion of sheared graphene sheets is studied with an in-house MD integrator in ref. 30. Energy dissipation during high speed sliding of graphene sheets is investigated with GROMACS^{31} and the DREIDING force field^{32} in ref. 33.
In addition to the relatively simple cases of nested nanotubes and graphene–graphene sliding, MD studies have evaluated more complex nanomachines. One example is a study of meshing gears made from functionalized carbon nanotubes^{34} using the Brenner potential.^{23} In ref. 18 a complex molecular planetary gear mechanism is analyzed using the UFF forcefield and Gasteiger partial charges. In ref. 35 a custom integrator is used with a custom force field based on CHARMM and AMBER force fields, with partial charges determined by ab initio calculations, to model the behavior of the experimentally realized “nanocars”^{3} rolling on a gold surface.
An experimental and MD study measured the energy required to force rotation of a sterically-congested molecular bond in a surface-bound molecule.^{36} The work to rotate one bond was 5 × 10^{−20} J. The energy barrier involved in rotating the bond in this molecule is much larger than that of the rotary joint presented in this paper. All of this energy is not necessarily dissipated since much could be recovered by a properly designed mechanism. Nevertheless, this study provides an upper bound on the energy dissipated, and also indicates the work that would need to be supplied and, ideally, recovered for a single rotation.
We model the effect of these torques with a phenomenological Langevin equation^{37,38} for the rotor's angular momentum L:^{2}
(1) |
The first term on the right-hand side is a frictional damping torque, with characteristic damping time 1/γ. The second term describes conservative torques, arising from a potential energy V(θ) depending on the rotor's orientation angle θ. For the machines considered in this paper, torques arising from the potential are relatively small. In the final term, W is a standard Wiener stochastic process^{37} modeling the fluctuating torques. Operationally, σdW denotes independent normally-distributed random values with zero mean and standard deviation over a time interval δt. The rotor's angular position changes as dθ/dt = ω, where ω = L/J is the angular velocity and J is the rotor's moment of inertia. Eqn (1) describes a stochastic process^{37} and numerical integration^{37,39} allows sampling solutions to this equation. This stochastic process is the rotational analog of a similar model applied to friction associated with linear motion.^{22,40}
The damping term in eqn (1) corresponds to a drag torque τ = γL = γJω. This drag is commonly expressed in terms of the rotational frictional drag coefficient, k_{rd}, defined by the relation τ ≡ k_{rd}ω, so that k_{rd} = γJ. This drag dissipates energy at the rate P_{dissipated} = k_{rd}ω^{2}. The energy dissipated due to this friction when rotating by angle ϕ at a uniform speed is τϕ. If this rotation takes place over time t, then ω = ϕ/t and the dissipation is
E_{dissipated} = P_{dissipated}t = k_{rd}ωϕ = k_{rd}ϕ^{2}/t | (2) |
When operating at temperature T, the fluctuation–dissipation theorem^{2,37,38} relates the damping and fluctuation parameters, γ and σ. The physical mechanism underlying this relation is that damping arises from continual interactions with the environment. Random variations in these interactions provide fluctuating torques that, on average, damp the rotor's motion. However, instead of coming to rest, these torques make the rotor continue to move with a fluctuating angular momentum that is zero on average but has a nonzero variance denoted by L^{2}_{th}, where L_{th} is called the thermal angular momentum. The stochastic process of eqn (1) has equilibrium variance σ^{2}/(2γ).^{37} Thus L^{2}_{th} = σ^{2}/(2γ). The thermal angular momentum is related to temperature through the equipartition theorem,^{41} which states the average rotational kinetic energy equals k_{B}T/2, where k_{B} is Boltzmann's constant. Since the rotational kinetic energy of a rotor with angular momentum L is L^{2}/(2J), the average kinetic energy in equilibrium is L^{2}_{th}/(2J). Thus L^{2}_{th}/(2J) = k_{B}T/2, so . Combining the expressions for L_{th} from the variance of eqn (1) and equipartition gives
σ^{2} = 2k_{B}TJγ | (3) |
In general, the drag can have a nonlinear dependence on the speed. Such nonlinearity can be significant at high speeds, e.g., close to the speed of sound in the structure or involving forces large enough to break bonds. On the other hand, at sufficiently low speeds, the linear term dominates in a Taylor expansion of the drag as a function of speed. As one comparison, the speed of sound in diamond is about 10^{4} m s^{−1}, whereas the examples in this paper involve rotation of carbon molecules with radius r ≈ 1 nm at up to angular velocities of ω = 100 rad ns^{−1}, with corresponding speed rω ≈ 10^{2} m s^{−1}.
The second assumption, uncorrelated random torques from other degrees of freedom, fails at sufficiently short time scales. For example, brief tilts of the rotor axis or oscillating normal modes of the housing could introduce short correlations in torques applied to the rotor. Eqn (1) assumes that these correlation times are short compared to times relevant for machine operation, which include the damping time. This requires a significant separation in time scales between the rotary degree of freedom and other motions. Molecular machines can achieve this separation with strongly bonded atoms throughout the structure except the parts intended to rotate, which instead have relatively weak interactions.
Thus we expect the stochastic model to apply to rotors operating at sufficiently low speeds and for sufficiently long times. We identify an appropriate time scale by comparing the standard deviation of changes in the rotor's angular momentum over various time scales δt. Specifically, consider ΔL(t) = L(t + δt) – L(t), the change in angular momentum starting at a time t. From a sample of the angular momenta, we compute s, the standard deviation of the changes ΔL at various times t in this sample, as a function of δt. Provided δt is short compared to the damping time 1/γ and the potential has only a small effect, the fluctuation term of eqn (1) dominates the changes, so . In this case, the ratio from the simulations should be nearly independent of δt. We evaluate this behavior with a “fluctuation plot” of this ratio vs. δt. The range of times δt achieving this independence indicates the time scale at which the molecular machine behaves approximately as a thermally-driven rotor.
Another check on the model validity is comparing the power spectrum of the angular momentum with that of the stochastic process of eqn (1) using parameters estimated from simulations of the rotor behavior. Frequencies with a close match indicate corresponding times over which the model is a good approximation. In our cases, we find the time scales determined from power spectra are similar to those from fluctuation plots.
The warmup stage starts by adjusting the initial geometry to an energy minimum without constraints on any part of the structure. Next, we impose boundary conditions, i.e., fix the positions of a specific subset of the atoms. This simulation starts the system at zero temperature and couples it to a temperature bath at 300 K. We run the warmup simulation long enough bring the atoms to this temperature.
The spin-down stage starts with a warmed up configuration and adds velocities to the atoms of the rotary component to give it a specified initial angular velocity ω_{0}. We then simulate the system in isolation, i.e., without coupling to a thermal bath. This is appropriate for the machines considered here, operating in vacuum, because the system is not in intimate contact with a thermal reservoir.
We repeat this procedure, with different random seeds, to obtain multiple samples of the behavior. Appendix A describes the simulation procedure in more detail.
On the other hand, fluctuations in the angular momentum are readily apparent in simulations over relatively short times. Thus we can use short simulations to estimate the fluctuation parameter σ and then determine the damping γ viaeqn (3). We do so by finding the value of σ that maximizes the likelihood of the observed changes in angular momentum over a fixed increment time δt according to eqn (1) and subject to the constraint of eqn (3). The second derivative of the log-likelihood at the maximum gives approximate confidence intervals for the estimated parameters.
In this fluctuation-based method, the choice of δt is somewhat arbitrary within a broad range. On the one hand, δt must be large enough that eqn (1) is a good approximation for the rotor, as discussed in section 3.2. On the other hand, δt must be short compared to the overall simulation time to have enough samples for estimation. As illustrated with the examples discussed below, this approach can estimate damping with simulation times short compared to the rotor's damping time even when operating at speeds comparable to thermal motion, so fluctuations are relatively large. For such situations, the fluctuation-based approach is computationally less costly than running simulations long enough to directly estimating damping with similar confidence.
In general, eqn (1) does not have a simple solution. In such cases, evaluating the likelihood requires sampling the stochastic process or numerically solving its Fokker–Planck equation.^{37} However, if the variation in the rotation potential V is small compared to thermal energies, the potential can be neglected, which greatly simplifies the estimation. In particular, when V is independent of the rotation angle, eqn (1) is an Ornstein–Uhlenbeck (O–U) process.^{42} In this case, the likelihood is readily computed. The examples considered in this paper have small potentials, allowing us to apply this simplification.
We applied the fluctuation-based method for evaluating friction to rotating nested nanotubes shown in Fig. 1. In this model system, a (9,9) single-wall carbon nanotube (SWCNT) rotates inside a (14,14) SWNCT. This example focuses on the rotary bearing itself, and ignores any damping that might arise from the housing to which the bearings would be connected in a molecular machine.
Fig. 1 Nested nanotubes. The inner tube is 6 nm long and the outer tube is 5 nm long. These consist of 882 and 1120 carbon atoms, respectively. |
The relatively small number of atoms in the nanotubes and their short damping time allow feasible simulations to show significant damping, even at relatively low rotation speeds, comparable to thermal speeds. This nanotube system thus provides a test case to compare damping estimated directly from the decrease in angular momentum with that estimated from fluctuations over a time scale much shorter than the damping time. Moreover, these simulations allow direct comparison with the extensive prior study of these nanotubes with simulations using other boundary conditions.^{21}Table 1 describes the inner nanotube, and Appendix A.1 describes our simulations.
Parameter | Value | |
---|---|---|
Temperature | T | 300 K |
Radius | r | 0.6 nm |
Mass | m | 1.76 × 10^{−23} kg |
Moment of inertia | J | 6.56 × 10^{−42} kg m^{2} |
Thermal angular momentum | L _{th} | 1.65 × 10^{−31} kg m^{2} s^{−1} |
Thermal angular velocity | ω _{th} | 25.1 rad ns^{−1} |
To identify the time scale over which the nanotubes match the stochastic model of eqn (1), Fig. 3 shows a fluctuation plot, described in section 3.2. The curves are roughly the same for δt above 0.01 ns. This indicates eqn (1) approximates the nanotube motion for times longer than this. In particular, short simulations, of just 0.1 ns, should be more than enough to estimate damping using the fluctuation-based method. These nested nanotubes have negligible potential barriers to rotation,^{26} so the potential is not relevant for friction estimation.
Method | Simulation time | Damping constant γ |
---|---|---|
Direct | 1 ns | 0.40 (±0.3)/ns |
Direct | 0.1 ns | 1.5 (±1.1)/ns |
Fluctuation | 0.1 ns | 0.45 (±0.2)/ns |
This estimate for nanotube drag is consistent with that of a previous study^{21} based on simulations at higher speeds, above 25 GHz. As described in Appendix A.1, this study gives k_{rd} = 3.4 (±0.8) × 10^{−33} kg m^{2} s^{−1}. This consistency indicates the linear fit of frictional force vs. speed of ref. 21 extends to the lower speed used in our simulation, i.e., 8 GHz.
We find similar behaviors when operating the nanotubes at higher speed and lower temperature, as described in Appendix A.2.
This section describes and evaluates such a design. The rotary joints consist of covalent acetylenic bonds between rotor and housing. Their linear geometry makes acetylenic bonds ideal as low barrier linkages in rigid organic frameworks. This property has been exploited in nanoscale mechanical design, where several designs of the Nanocar^{3} employed acetylenic linkages as axles between fullerene wheels and a planar chassis. The structure includes enough of a housing around the rotor to include long range interactions between the rotor and the rest of a molecular machine using the rotary joint. As with the nanotube bearing, we consider operation in vacuum at room temperature.
To evaluate rotor behavior, we used two computational techniques. First, we applied density functional theory (DFT) calculations to the rotor and part of the surrounding structure under symmetry constraints to determine how the rotor interacts with the housing. This included the predicted rotational barrier due to the shortest range rotor-housing interactions, and the vibrational mode energies associated with the rotor motion. Second, we used the DFT results to generate missing AMBER force field terms and RESP partial charges of the rotor and surrounding structure for the molecular dynamics simulations, which then were used to sample the rotor's motion within the full structure. We used these samples to estimate the drag. Appendix A.3 describes our simulations.
Parameter | Value | |
---|---|---|
Temperature | T | 300 K |
Mass | m | 4.24 × 10^{−24} kg |
Moment of inertia | J | 1.76 × 10^{−42} kg m^{2} |
Amplitude of potential | V _{0} | 4.15 × 10^{−22} J |
Thermal angular momentum | L _{th} | 8.54 × 10^{−32} kg m^{2} s^{−1} |
Thermal angular velocity | ω _{th} | 48.5 rad ns^{−1} |
We designed this specific molecule as likely to give a rotary joint with low dissipation. Specifically, this molecule:
• provides a complete rotor and housing that includes all thermal motion in and around the rotor, and long-range interactions between the rotor and its housing
• isolates the rotor motion from the housing vibrations by use of a stiff material (diamond)
• has a low rotational energy barrier
• has closely matched rotor length and housing spacing so there is little strain (compression or stretching) between the C_{2} spacers and either the rotor or the housing.
This molecule has two properties that simplify computational studies. First, the high symmetry of the rotor and nearby parts of the housing reduces the computational cost of determining its properties, including geometries corresponding to extrema of the rotational potential, normal modes of vibration, and partial charges on the atoms. Second, the molecule consists of just three atom types: sp3 (tetrahedral) carbon, terminal C–H bond hydrogen, and sp (acetylenic) carbon. These standard atom types readily transfer to any molecular force field for additional testing and validation.
We use these charges as part of the force field for our molecular dynamics simulations, both to sample rotor behavior for determining drag and to estimate the rotation potential V(θ). This potential is an aggregate property of the molecular dynamics force field. It quantifies the conservative torques on the rotor in the stochastic model given by eqn (1). As described in Appendix A.5, the potential energy as a function of rotation angle is well-approximated by
V(θ) = V_{0}(1 − cos(3θ)) | (4) |
This potential barrier, 2V_{0} = 8.3 × 10^{−22} J, is comparable to estimates of barriers for rotation about triple bonds.^{2} This barrier is about 0.2k_{B}T, which is relatively small compared to thermal energy.
The potential applies a torque –dV/dθ to the rotor. The maximum torque is 3V_{0} = 1.2 pN nm. As an estimate of corresponding forces, if the maximum torque is mostly at the rotor's outer rim, with radius r = 0.9 nm (see Fig. 11 in Appendix B), the force is 3V_{0}/r = 1.4 pN. If instead the torque mainly acts through the axle, of radius r = 0.05 nm, the force is 25 pN.
The atoms of the rotor change position slightly during the simulations. However, due to the strong bonding, these motions have only a small effect on the moment of inertia. In particular, the standard deviation of the moment of inertia about the rotor axis, J, is 6 × 10^{−45} kg m^{2}, less than 1% of the mean value, given in Table 3.
To identify the time scale over which the rotor matches the stochastic model of eqn (1), Fig. 6 shows a fluctuation plot, described in section 3.2. The values are roughly the same for δt above 0.2 ns or so. For shorter times, the observed fluctuations are larger than the expected behavior from eqn (1). This indicates the stochastic model does not correctly describe short-time motion of the rotor, below about 0.2 ns or above 5 GHz. This conclusion is similar to that from the peaks in the power spectrum shown in Fig. 7. These observations indicate eqn (1) provides a reasonable approximation of the rotor's motion for times longer than about 0.2 ns.
Fig. 6 Standard deviation of changes in rotor angular momentum in time δt divided by , for each of the rotor simulations (upper, colored curves). The dashed and solid black curves show the value of this ratio expected from the stochastic process model of the rotor, eqn (1) discussed in section 5.5, with and without the potential, respectively. Times are shown on a logarithmic scale. |
Fig. 7 Average of the angular momentum power spectra from the spin-down simulations. The horizontal axis shows frequency on a logarithmic scale. Power values are in decibels, normalized so the total power equals the sum of the squared angular momenta from the simulations. For comparison, the dashed and solid black curves show the power spectrum of the stochastic model of the rotor, eqn (1), with and without the potential, respectively. |
The potential is relatively small compared to k_{B}T. In addition, δt above 0.2 ns involves motion of more than a complete revolution of the rotor for the angular velocities in the simulation, which remain close to the initial 100 rad ns^{−1} because damping does not slow the rotor much during the 5 ns simulation time. Thus the small effects of the potential are further reduced by the rotor averaging over the potential during the time steps we consider.
For example, with the angular velocity in our simulations, ω_{0} = 100 rad s^{−1}, the time for a full rotation is t = 2π/ω_{0} = 0.06 ns. From eqn (2), the frictional energy dissipation during that rotation is 2πk_{rd}ω_{0} ≈ 1.5 × 10^{−23} J. By comparison, a full rotation of this molecular rotor moves through three potential wells, each of depth 2V_{0} (see Table 3). Thus climbing the potential barriers requires 6V_{0} ≈ 2.5 × 10^{−21} J. However, this work could, in principle, be recovered as the rotor moves into each potential well. The energy dissipated by friction is about 100 times smaller than the (possibly recoverable) energy required to climb the potential, which in turn is an order of magnitude smaller than the experimentally measured energy, 5 × 10^{−20} J, to force rotation of a surface-bound molecule.^{36}
Fig. 6 shows one such comparison. The rotational potential accounts for some of the increase in fluctuations at short times. In particular, the rotor converts between kinetic and potential energy as it moves between the top and bottom of the potential wells. For our simulations, the angular velocity remains close to its initial value 100 rad ns^{−1}. The potential, in eqn (4), is 3-fold symmetric, so rotation from one potential well to the next requires rotation through 2π/3 ≈ 2 rad, which occurs in 0.02 ns, corresponding to 50 GHz. This motion leads to the oscillating fluctuation size shown in the dashed curve of Fig. 6.
Another perspective on short-time rotor behavior is the average of the simulations' power spectra, shown in Fig. 7. By comparison, the power spectrum of an O–U random process is proportional to 1/(γ^{2} + (2πf)^{2}).^{37} For frequencies large compared with γ, this power decreases as 1/f,^{2} and hence appears as a straight line on a log–log plot. Fig. 7 compares this behavior with that of the rotor simulation, showing close correspondence with the rotor's angular momentum up to about 10 GHz.
The behavior including the potential (dashed curve in Fig. 7) has the broad peak near 50 GHz, showing the potential ripple accounts for both the frequency (due to time it takes rotor to move over one potential barrier at 100 rad ns^{−1}) and the magnitude of the peak.
Above about 50 GHz, both the fluctuation plot and power spectrum show deviations from the stochastic process model, even when including the potential. One source of these deviations is temporal correlations in torques applied to the rotor, e.g., due to movement of its axis of rotation. The oscillation frequencies of the structure's normal modes indicate the time scale of these motions. We evaluated the normal modes of the rotor-housing structure using GROMACS. The lowest mode is the rotor twisting about its axis, at about 7.24 GHz. This is close to the oscillation frequency of the potential minimum: . This mode corresponds to the rotary motion given by the stochastic model. The next lowest modes involve the rotor tilting, at around 100 GHz. These modes characterize deviations from the stochastic model, and show the large separation in time scales between the rotary degree of freedom and other modes in the structure. This separation is the expected behavior for a rotary bearing consisting of strongly bonded structures and weak barriers to rotation.
These comparisons between simulations and the stochastic process with the estimated parameters show the atoms in the rotor and housing behave approximately as a stochastic rigid rotor, with a single degree of freedom around the rotary joint. The atoms' other degrees of freedom act as a source of uncorrelated random fluctuations, i.e., a thermal bath.
The method was applied to nested nanotubes and a rotor covalently bonded to a housing by single acetylenic bonds. For the nested nanotubes, damping time is 1/γ = 2.4 ns and drag coefficient is k_{rd} = 2.9 × 10^{−33} kg m^{2} s^{−1} (section 4.2). The bonded rotor joint has 1/γ = 70 ns and k_{rd} = 2.4 × 10^{−35} kg m^{2} s^{−1} (section 5.4). From eqn (2), this means rotational friction dissipates 120 times less energy for the bonded rotor than for the nested nanotubes when they rotate at the same speed. Thus acetylenic rotary joints can have exceptionally low dissipation, which could make them useful bearings in future molecular machines.
Machine design often focuses on performance and ease of manufacture. For molecular machines that are not yet possible to fabricate, an important additional criterion is the computational cost to simulate their behavior. As an example of designs that simplify drag analysis, rotors with small potential barriers compared to thermal energy allow neglecting the potential in the stochastic model used with the fluctuation-based method. In addition, strong bonding among degrees of freedom other than those of the machine's operation (in this case the single rotational degree of freedom) allows treating the machine as a rigid body, with other degrees of freedom providing uncorrelated random fluctuations.
An open question is the range of conditions under which the fluctuation-based method performs well. For instance, potentials with larger barriers to rotation, less rigid molecular housings or less symmetric rotational geometry than the examples considered in this paper may significantly reduce the accuracy of modeling the rotor as a stochastic process with a single degree of freedom. Since neglecting the rotational potential simplifies the analysis, a specific issue for future work is determining design rules under which a molecular rotor is sufficiently well-approximated as having zero potential when evaluating its friction.
The fluctuation-based method discussed in this paper could evaluate drag in a range of operating conditions. Examples are how rotary friction depends on rotor size, rotation speed and temperature. This could require accounting for additional physical effects. For instance, when operating well below room temperature, the variation in the potential is significant compared to k_{B}T, leading to more complicated behavior than seen in the examples of this paper.^{2,43} Moreover, at low temperatures, quantum effects alter how energy spreads among degrees of freedom and introduce additional dissipative mechanisms, as described in Appendix A.6.
The interactions between rotor and housing arise mainly from forces between the rotor and nearby atoms in the pyramid connectors. These interactions largely determine the potential barrier to rotation. Future work could more precisely quantify these short-range interactions and the extent to which they are responsible for the friction. Identifying the main contributions to friction will suggest how changes to the housing structure affect dissipation, e.g., by altering stresses on bonds near the rotary joint.
The fluctuation-based method can apply to many molecular machines in addition to the rotary bearings considered in this paper. In particular, a single bearing will only be a small part of a useful molecular machine using the bearings, such as a computer.^{20} Designing such machines, containing large numbers of atoms, will require efficient computational methods. Thus the fluctuation-based method could aid computational explorations of a wide variety of designs to identify low-power molecular machines. Fluctuations can estimate drag from simulations much shorter than damping time, which is particularly useful for low-friction molecular machines with many atoms. However, the simulation must still be long enough to average over dynamics of other modes of behavior than the degree(s) of freedom involved in the machine's intended operation. Otherwise fluctuations may not be independent, as assumed by the stochastic process model. More specifically, the method is useful when the correlation time of fluctuations is much smaller than the damping time of the machine's operation. Molecular machines made of strongly bonded structures can have such separation of time scales, as illustrated by the examples discussed in this paper. Thus the method described in this paper could be a useful tool for evaluating dissipation in such machines.
For both the nanotube and rotor simulations, the warmup stage starts with all atoms at zero velocity and couples them to a heat bath, i.e., the simulations used a canonical (NVT) ensemble. The spin-down simulations used a microcanonical (NVE) ensemble, i.e., without temperature control.
The warmup and spin-down simulations are performed in vacuum and use 1 fs time steps. In accordance with best-practices for molecular dynamics simulation in vacuum while reducing computational resource use, all GROMACS simulations were performed under periodic boundary conditions using 25 × 25 × 25 nm boxes with rlist, rcoulomb, and rvdw cut-offs of 10 nm, particle mesh Ewald (PME)^{45} coulomb type and vdwtype (with pme-order of 6), and the Verlet cutoff scheme under a grid-based neighbor list determination.
The box size used for periodic boundary conditions is larger than the sum of half the molecule's dimensions and the cut-off lengths. The volume of this box remains constant during the simulation. Since the molecule is in vacuum, the simulation has both constant volume for the full simulation volume and zero pressure on the molecule.
The warmup simulation uses the following parameters: 50000 steps at 1 fs per step (50 ps total simulation time), xyz periodic boundary conditions with Verlet cutoff, 10 nm cutoff distance (which exceeds the largest atom-to-atom distance in the structure, and is less than the distance between copies of the structure due to periodic boundary conditions), v-rescale temperature control with 1 ps time constant, and random seed set to 42, 43, 44, 45, or 46. These seeds give five different samples of the structure at 300 K. During the warmup simulations, the structures reach 300 K within about 10 ps, indicating the 50 ps warmup simulation time is more than sufficient to eliminate transients during warmup.
Table 5 provides the force field parameters^{46} used for the nested nanotube bearing simulations. The bond stretch and bond angle terms are based on the AMBER force field.^{47} The Ryckaert–Bellemans dihedral terms are based on the bending rigidity of graphene,^{48} and the Lennard-Jones terms are based on ref. 27.
Interaction type | Parameter | Value | Unit |
---|---|---|---|
Morse potential bond stretch | b | 0.142 | nm |
Morse potential bond stretch | D | 722.7 | kJ mol^{−1} |
Morse potential bond stretch | β | 18.33 | nm^{−1} |
Harmonic angle potential | k | 527.2 | kJ mol^{−1} rad^{−2} |
Harmonic angle potential | θ _{0} | 120.0 | Degrees |
Ryckaert–Bellemans dihedral | C_{1}, C_{3}, C_{4}, C_{5} | 0.0 | kJ mol^{−1} |
Ryckaert–Bellemans dihedral | C _{0} | 17.843 | kJ mol^{−1} |
Ryckaert–Bellemans dihedral | C _{2} | −17.843 | kJ mol^{−1} |
Lennard-Jones | ∈ | 0.28454 | kJ mol^{−1} |
Lennard-Jones | σ | 0.341 | nm |
Our nanotube simulations are similar to the “coast-down” method of ref. 21. Minor differences arise from our use of GROMACS instead of LAMMPS, the force field of Table 5 instead of AIREBO, and having the outer tube anchored at one edge instead of freely counter-rotating.
Our friction estimate is consistent with the linear fit to normalized values given in ref. 21. This fit is F* = Av* where F* and v* are normalized friction force and sliding velocities [eqn (23) and (24)],^{21}A = 4.62 × 10^{−4} and the fit is consistent with having zero friction at zero velocity.† Converting the friction force to a torque, and the sliding velocity to the corresponding angular velocity, the linear relation between F* and v* gives the rotational frictional drag coefficient described in section 3.1:
(5) |
When the nanotube rotates at speeds well above thermal velocity, fluctuations are relatively unimportant. In that case, the decrease in angular momentum is readily apparent, as illustrated in Fig. 8 for nested nanotubes with initial angular velocity 500 rad ns^{−1}. This initial speed is ten times the angular velocity considered in section 4, and is well above typical thermal speeds (see Table 1). Table 6 compares drag estimated with the direct and fluctuation-based methods for this case. Because thermal fluctuations are relatively small, the direct method gives good estimates even over short times. Both methods give comparable drag estimates, and are consistent with the values at the lower speed, shown in Table 2. This similarity is consistent with the linear relation between frictional torque and angular velocity found in other studies,^{21,27} as well as the linear form for the drag in eqn (1).
Method | Simulation time | Damping constant γ |
---|---|---|
Direct | 0.5 ns | 0.38 (±0.03)/ns |
Direct | 0.1 ns | 0.46 (±0.1)/ns |
Fluctuation | 0.1 ns | 0.41 (±0.1)/ns |
Fluctuations and damping depend on temperature. As an example, Table 7 shows the damping estimates from the two methods at T = 80 K. We performed 5 simulations for each of the two initial angular velocities shown in the table. In this case, the thermal angular velocity is ω_{th} = 13.0 rad ns^{−1} so even the smaller initial speed, ω_{0} = 50 rad ns^{−1}, is well above the thermal speed. This means fluctuations are less important than at higher temperatures, and both the direct and fluctuation-based methods give similar results, as shown in Table 7. At 80 K, the damping time 1/γ is about 15 ns, so there is relatively little damping during our 0.5 ns simulations.
Method | Initial angular velocity | Damping constant γ |
---|---|---|
Direct | 500 rad ns^{−1} | 0.067 (±0.01)/ns |
Fluctuation | 500 rad ns^{−1} | 0.066 (±0.01)/ns |
Direct | 50 rad ns^{−1} | 0.071 (±0.06)/ns |
Fluctuation | 50 rad ns^{−1} | 0.068 (±0.03)/ns |
The fluctuation–dissipation theorem, eqn (3), relates the fluctuation (σ) and damping (γ) parameters to temperature (T). Specifically the ratio σ^{2}/γ is proportional to T. However, this does not determine how the two parameters vary individually. Comparing the damping at 80 K and 300 K identifies how these parameters depend on temperature. Fitting the parameter estimates from our simulations at both temperatures to a power-law dependence, σ ∝ T^{α} and γ ∝ T^{β}, gives α = 1.23 ± 0.08 and β = 1.47 ± 0.15. This is consistent with a study of nanotube friction over a wider range of temperatures that found β = 1.53 ± 0.04.^{27}
Structure optimizations, normal mode analyses, and molecular dynamics simulations used the AMBER99 (ref. 47) force field C–H and C–C bend, stretch, dihedral, and C6/C12 parameters serving as a basis for the calculations. Terms for the CC (triple bond) and C–C# bonds were defined as given in ref. 50.
The full structure (rotor plus housing) is first energy-minimized without applying anchors to any part of the structure. Next, anchors are applied to 514 atoms in the stem of the housing, and the structure is heated to 300 K in a warmup simulation.
The output of each of the five warmup simulations is a molecular structure in G96 format which contains the position and velocity for each atom in the structure. The initial angular velocity of the rotor is set by using a Matlab/GNU Octave script to adjust the velocities of individual atoms. First the principal axes of the rotor are determined, and then the full system is rotated into a coordinate frame in which the z axis aligns with the rotor's main principal axis. The component of the initial angular momentum due to thermal excitation about the main principal axis L^{warmup}_{z} and initial polar moment of inertia about the main principal axis J_{z} are calculated. A correction velocity of
(6) |
The spin-down simulation uses the speed-adjusted G96 files as inputs and anchors the atoms in the stem. The simulation parameters are: 5000000 steps at 1 fs per step (5 ns total simulation time), xyz periodic boundary conditions with Verlet cutoff, 10 nm cutoff distance.
Table 8 shows the energy drift for the 5 ns spin-down simulation with random seed = 42, calculated with the GROMACS utility g_energy_d.^{46} The other simulations have similar values. For the N = 9831 atoms free to move, 3Nk_{B}T = 1.2 × 10^{−16} J (or 7.4 × 10^{4} kJ mol^{−1}), about half the average energy in the table. This is due to GROMACS including an offset of about 9.46 × 10^{4} kJ mol^{−1} from the minimum-potential configuration of the rotor.
Quantity | Average | Standard deviation | Total drift | Units |
---|---|---|---|---|
Total energy | 1.7 × 10^{5} | 12 | 3.3 | kJ mol^{−1} |
Temperature | 303 | 1.8 | −0.01 | K |
During each of the five spin-down simulations, the total system energy remained nearly constant. Specifically, over the full five nanosecond spin-down, the largest increase in total system energy, as computed by the GROMACS analysis tool g_energy_d, of any of the simulations was 6.51 × 10^{−21} J (4 kJ mol^{−1}). This drift is small enough to not effect the conclusions of this study.
Fig. 9 Potential energy of the rotor and housing structure as a function of rotor rotation angle within one of the three symmetric potential wells. The points are the values measured with GROMACS and the curve is eqn (4). The angle 0° is the configuration shown in Fig. 4. |
For the purposes of our analysis, the key quantity is the height of the potential barrier 2V_{0} compared to k_{B}T. As a check on the size of the potential, we performed a spin test at zero temperature. This test started the rotor spinning and measured the change in angular momentum as the rotor moved from one potential well to the next. This case has no thermal motion so the rotor is affected only by the potential. The observed changes were consistent with eqn (4) and the value V_{0} in Table 3.
Fig. 10 shows the configuration of atoms at the rotary joint corresponding to the extremes of the potential. The staggered geometry is the minimum of the rotational potential, while the eclipsed geometry is the maximum. The staggered configuration corresponds to rotor rotation angles (see Fig. 4) of θ = 0, 120, 240 degrees. The eclipsed configuration corresponds to angles of θ = 60, 180, 300 degrees.
Fig. 10 Ball-and-stick and space-filling models of the molecular joint in two configurations: eclipsed and staggered on left and right, respectively. C_{2} spacers are indicated by the arrows. These geometries have molecular D_{3d} and D_{3h} symmetries,^{52} respectively. |
The potential arises mainly from non-bonding interactions between atoms of the rotor and the pyramid portions of the housing. The single C_{2} spacer is, then, a compromise between minimizing the rotational barrier and maintaining stiffness along the rotation axis.
The thermal de Broglie wavelength of an object with mass m is which is 2 × 10^{−3} nm for the rotor, much smaller than its size. Thus the rotor is well-approximated as a localized classical mass.
At the temperature considered here, typical angular momenta of the nanotube and rotor (see Tables 1 and 3) are about 1000 times larger than the quantum of angular momentum ℏ = 10^{−34} kg m^{2} s^{−1}. Thus the classical approximation of continuous variation of angular momentum is reasonable. These observations indicate the rotor behaves as a classical object for operation at up to at least a few GHz.
Drag arises from spreading organized energy in rotor motion to the many other degrees of freedom of the molecular structure. Thus in addition to classical physics being sufficient to describe the rotor motion itself, molecular dynamics simulations assume classical behavior for other modes of motion. In particular, classical equipartition of energy among these modes only applies when thermal energy, k_{B}T, is large compared to the quantum level spacing hf, where f is a mode's frequency.
The rotor speeds we consider (below 100 GHz) have k_{B}T large compared to the quantum level spacing. Thus the classical behavior assumed in molecular dynamics simulations is a reasonable approximation for the rotor itself.
In addition to the rotary motion, the simulations include much higher frequency modes, particularly the hydrogen bond motions. Resolving these motions requires simulation time steps, 1 fs, much shorter than the time scales of the rotary behavior. For these high-speed motions, the quantum energy spacing is large compared to k_{B}T. Thus these modes will have less energy than predicted by classical equipartition and GROMACS simulations. As a check on whether this suppression affects the rotor drag results, we repeated the simulations with the GROMAS H-bond constraint, which effectively keeps the hydrogen–carbon bonds in their ground state throughout the simulation. With this constraint, we find the same value for the rotor drag, within the reported confidence intervals, as in our original simulations. Nevertheless, this constraint significantly reduced the short-time fluctuations. With this reduction, the fluctuations closely match those expected from the stochastic process model (shown as the dashed curve in Fig. 6) for time differences as short as 1 ps. This correspondence suggests the large fluctuations at short times seen in Fig. 6 may be a simulation artifact, arising from ignoring quantum suppression of high frequency modes.
Quantum effects are more important at low temperature, where thermal energy is smaller. Furthermore, dissipative forces arising from quantum effects^{53} may be important at low temperatures.^{54}
The rotor itself is a hexagonally-shaped piece of hexagonal-lattice diamond, or lonsdaleite, as shown in Fig. 11. DFT calculations indicate that lonsdaleite would be stable,^{58} and it has been fabricated.^{59} As with the housing, the hydrogen-passivated facets of the rotor are predicted to be stable and free from surface reconstruction.
A ball-and-stick atomic model of the rotor and its two rotary joints is shown in Fig. 13. Each rotary joint is composed of four co-linear carbon atoms, similar to the organic molecule 2-butyne. The single bonds at the outside ends of the co-linear chain (e.g., bond C9442–C9443 and bond C9444–C9446) provide for the rotation. Such four-carbon chains have been shown experimentally to function as stable mechanical rotary joints.^{3} Of the four carbon atoms, the outside two (e.g., C9422 and C9446) are embedded within the framework of the rotating components – one within the housing and one in the rotor, while the inner two (e.g., C9443 and C9444) are free to rotate. We refer to the inner two carbon atoms as a “C_{2} spacer”.
Many proposed applications for molecular machines based on the rotary joint structure are intended to operate in vacuum (e.g. mechanical computation devices). Nevertheless, it is likely that the structure will be stable at room temperature in atmosphere. Several examples of atomically precise mechanical interactions are known to be stable in air, including the self-retraction of nested nanotubes,^{14} self-retraction of graphene flakes,^{60} and rotation of acetylenic axles.^{3} In addition, the hydrogen-passivated carbon {111} surface is chemically stable in air.^{61,62} These experimental results on similar systems strongly suggest that the rotary joint molecule would be stable under ambient conditions.
Fig. 12 Rotor-pyramid assembly used for DFT calculations, shown as insets and in the context of the full structure. |
The rotor and housing structure conforms to C_{2v} symmetry^{52} at certain rotor geometries. In the absence of the housing, the rotor-pyramid assembly has significantly higher symmetry at certain rotor geometries. Specifically, defined against the nearest-neighbor H⋯H interactions between the rotor and pyramids, two symmetry-constrained geometries bring the 12 total pairs of H⋯H interactions between rotor and pyramid eclipsed and staggered conformations with D_{3d} or D_{3h} symmetries, respectively. Fig. 10 shows these conformations. These symmetries allow calculations at significantly higher (i.e., more accurate) levels of theory than is possible for lower-symmetry structures.
The B3LYP hybrid density functional^{64,65} was employed for all calculations, along with various calculations with the 6-31G(d,p),^{66} 6-311G(2d,p), and 6-311G(2d,2p)^{67} Gaussian-type basis sets and D3-type Grimme Dispersion correction.^{68}
RESP model^{69,70} charge calculations for the rotor-pyramid assemblies were produced at the B3LYP/6-31G(d,p) level of theory, employing the Merz–Singh–Koller scheme^{71,72} as a basis for Antechamber script input (program option “Pop=MK IOp(6/33=2,6/41=10,6/42=17)”).
To reduce the computational cost of the DFT calculations, we restricted the evaluation to rotor-pyramid assemblies shown in Fig. 14. These assemblies have few enough atoms to allow feasible DFT calculations, while including most of the interactions between the rotor and the housing. We further reduced the cost by only evaluating the highly symmetric eclipsed and staggered conformations of the assemblies described in Appendix C.1. These conformations correspond to the extrema of the potential evaluated with GROMACS in Appendix A.5. Evaluating just these two cases is sufficient to compute the potential barrier, i.e., as the difference in energy between the eclipsed and staggered geometries of the joint.
Table 9 shows the values for the barrier, i.e., 2V_{0}, for assemblies with various numbers of atoms in the pyramid. The calculated barrier between eclipsed and staggered geometries is well below k_{B}T = 2.5 kJ mol^{−1} at the temperature we consider. The rotational barrier does not change noticeably when including more atoms from the pyramid in the calculation, indicating that the barrier is a local phenomenon primarily due to interactions among atoms close to the rotary joint. The small size of these computed energy barriers (0.05 kJ mol^{−1}) are in excellent agreement with experimental measurements and theoretical calculations for the rotational barrier around the acetylene linkage in the symmetric 2-butyne molecule.^{74} Diphenylacetylene, a structurally similar molecule with closer H⋯H interactions than diadamantylacetylene and a preference for maintaining planarity, has a measured rotational barrier of 2.5 kJ mol^{−1} (calculated to be 3.3 kJ mol^{−1} at the B3LYP/6-311+G(d,p) level of theory [Table 2]).^{74} At the operational temperature we consider, the barrier for the C_{2} spacer motif is expected to be inconsequential to the operation of the full assembly.
PyrN-C_{2}–Rot–C_{2}-PyrN | ||||
---|---|---|---|---|
DFT theory level | C_{10}–C_{2}–C_{10} | Pyr1 | Pyr2 | Pyr3 |
B3LYP-GD3/6-311G(2d,2p) | 0.022 | — | — | — |
B3LYP-GD3/6-31G(d,p) | 0.030 | 0.063 | −0.006 | 0.096 |
B3LYP/6-31G(d,p) | 0.003 | 0.122 | 0.057 | 0.154 |
The barrier computed by GROMACS for the complete structure, 2V_{0} = 0.5 kJ mol^{−1} from Table 3, is larger than that determined by DFT for isolated rotor-pyramid structures (−0.006 to +0.154 kJ mol^{−1}, see Table 9). Molecular mechanics force fields using RESP charges are of inherently lower accuracy than DFT methods.^{47} It is therefore not surprising that GROMACS overestimates the rotational barrier since the barrier predicted by DFT is extremely small. The larger barrier calculated by GROMACS may be due to interactions in the full structure modeled by GROMACS that are not present in the smaller structures modeled by DFT, such as strain induced by the housing or long range interactions between RESP charges. Regardless of the origin of the discrepancy, a barrier height of 0.5 kJ mol^{−1} is comparable to values reported for smaller molecules with similar rotational structure, e.g., 12 kJ mol^{−1} for internal rotation in ethane and 0.07 kJ mol^{−1} for 2-butyne.^{1,74}
Footnote |
† The corrected linear relation between force and speed^{49} corresponds to [Fig. 9]^{21} rather than [eqn (28)].^{21} |
This journal is © The Royal Society of Chemistry 2017 |