Molecular dynamics simulation of thermal ripples in graphene with bond-order-informed harmonic constraints

Atomistic molecular mechanics, including molecular statics (MS) and molecular dynamics (MD), has become one of the most widely used particle-based simulation methods in nanoscale physics and chemistry. In particular, it has contributed valuable insight into the properties of carbon nanotubes and graphene, which have enjoyed significant attention in the last decade. Some of the important results from molecular mechanics include the prediction of exceptionally high thermal conductivity in $s{{p}^{2}}$-bound atomically thin carbon shells [1], a numerical evaluation of graphene's Young's modulus [2] and bending rigidity [3] and an understanding of phonon scattering at the graphene-fluid interface [4].

The interatomic interactions that describe local binding provide the core of the physical representation of simulated systems. One method of describing interatomic interactions for graphene and carbon nanotubes that is well established in terms of overall accuracy and relative computational efficiency is the empirical bond-order Brenner potential, which is specifically parameterized to reproduce a variety of experimental data, including the phonon spectrum of graphene [5–7]. Numerous works have also employed the highly parameterizable modified embedded-atom method (MEAM), which is formulated to include both the bond pair interactions and the local electron density contribution, mimicking the density functional theory approach [8–12]. These highly successful descriptions of interatomic interactions in covalent materials are relatively computationally intensive and have been implemented in many simulation packages.

A great variety of considerably less computationally intensive representations have existed for several decades. In particular, a number of widely used molecular dynamics simulation packages rely on force fields, which utilize locally decoupled harmonic geometry constraints in an atomistic representation of molecules [13, 14]. Such a representation is highly computationally efficient and was originally designed for small organic molecules. Despite the wide availability of this model, successful large-scale dynamic simulations of atomically thin layers with harmonic constraints have been relatively rare, partly due to their inability to describe bond formation and breaking and their lack of reasonable parameterization for intact membranes. A good starting point for determining suitable harmonic constraints for intact graphene membranes was provided in [15], where the bond and in-plane angular constraint stiffness values were analytically obtained from the bond-order potential for a structural model of graphene, followed by quasistatic tests on small rectangular samples. The dynamic performance of harmonic constraint-based interactions in considerably larger systems generally remains an open question.

The performance of the harmonic constraint model in describing the dynamic flexural waves in graphene is of particular interest. The heights of these ripples have been theoretically predicted to be comparable to [16–18] and were recently measured to exceed [19], the lattice constant of graphene. Moreover, the effective motion time scales for these ripples are well within the nanosecond to microsecond simulated times [18] that are currently routine for MD. Therefore, one can expect such dynamic flexural vibrations to have a considerable effect on simulated systems, e.g. on graphene's interaction with fluids [4] and gases [20].

In this work, we use a complete set of bond, in-plane angle and dihedral angle constraint constants directly obtained from a bond-order potential and test the resulting force field against the results obtained with the use of the bond-order potential as well as the MEAM. We do so by dynamically simulating and statistically quantifying thermally induced ripples in relatively large graphene sheets at various temperatures. We demonstrate that despite presenting an intrinsically stiffer model of graphene in terms of the overall amount of thermal rippling (as compared with the bond-order model and previously published models), the harmonic constraint model is, overall, capable of representing the dynamic morphology of graphene reasonably well. Our results are demonstrated for various temperatures and membrane sizes. We also observe that within the harmonic constraint model, the out-of-plane ripple relaxation occurs via the activation of long flexural modes.

The interatomic interactions described by the bond-order potential (including the first-generation [5] and the second-generation [7] Brenner potential) define the following many-body expression for the covalent bond energy between the nearest neighbors i and j:

$$\begin{eqnarray}&&{{V}_{ij}}={{V}_{R}}\left( {{r}_{ij}} \right)-{{B}_{ij}}\left( \mathop \sum \limits_{k\ne i,j} G\left( {{\theta }_{ijk}} \right),\mathop \sum \limits_{k,l\ne i,j} H({{\alpha }_{ijkl}}) \right)\;{{V}_{A}}\left( {{r}_{ij}} \right),\;\end{eqnarray} \tag{ 1 }$$

where ${{V}_{R}}\left( {{r}_{ij}} \right)$ and ${{V}_{A}}\left( {{r}_{ij}} \right)$ are the repulsive and attractive Morse-type [21] i-j pairwise distance terms, respectively, and ${{B}_{ij}}\left( \sum G,\;\sum H \right)$ is the bond order functional, which depends on the in-plane inter-bond angles ${{\theta }_{ijk}}$ formed by the i-j and j-k bonds and the dihedral angles ${{\alpha }_{ijkl}}$ between the planes formed by the i-j-k and i-j-l particle triplets. Equation (1) thus explicitly includes the local bonding environment of the i-j pair, which affects the bond lengths and angles in a mutually dependent manner.

The harmonic geometry constraint model, on the other hand, presents a considerable simplification compared with equation (1). The total potential energy of the simulated system, with the total number of bonds, in-plane angles, and dihedral angles of ${{N}_{bond}}$, ${{N}_{ang}}$, and ${{N}_{dih}}$, respectively, is:

$$\begin{eqnarray}\begin{array}{lllllllllllllll} {{E}_{tot}} & = & \frac{1}{2}{{k}_{bond}}\mathop \sum \limits_{{{N}_{bonds}}} {{\left( {{r}_{ij}}-a \right)}^{2}} \\ {} & {} & +\frac{1}{2}{{k}_{ang}}\mathop \sum \limits_{{{N}_{ang}}} {{\left( {\rm cos} {{\theta }_{ijk}}-{\rm cos} {{\theta }_{0}} \right)}^{2}} \\ {} & {} & +\;\frac{1}{2}{{k}_{dih}}\mathop \sum \limits_{{{N}_{dih}}} {{\left( {{\alpha }_{ijkl}}-{{\alpha }_{0}} \right)}^{2}}, \\ \end{array}\end{eqnarray} \tag{ 2 }$$

where $a$, $\;{{\theta }_{0}}$, and ${{\alpha }_{0}}$ are the prescribed bond length, in-plane angle and dihedral angle, respectively. In the case of graphene, $a=1.42\;{\rm {\ \mathring{\rm A} }}$, and $\;{{\theta }_{0}}=\frac{2\pi }{3}\;\;{\rm rad}$. A total of four dihedral energy contributions per bond are present in two pairs with ${{\alpha }_{0}}=0\;{\rm rad}$ and ${{\alpha }_{0}}=\pi \;{\rm rad}$, respectively. The stiffness constants ${{k}_{bond}}$, ${{k}_{ang}}$, and ${{k}_{dih}}$ are the only effective fitting parameters in equation (2)—a small number compared to the bond-order potential or the MEAM, which, depending on the form, can contain more than ten parameters. Granted, the presence of only three adjustable parameters, especially within the representation with limited bond-angle coupling, limits our ability to fine-tune equation (2) to fit a large set of empirically known material properties. However, as we show here, a relatively accurate representation of a thermally fluctuating graphene membrane is possible compared with the bond-order potential and the MEAM. Considering the low computational cost and high availability of harmonic constraints as a generic MD model, harmonic constraints may be an attractive choice for large-scale multicomponent simulations.

One can see that in order to mimic the bond-order representation, the stiffness constants in equation (2) can be directly obtained from the total potential energy of the system, represented by an arbitrary detailed interaction potential via evaluating the corresponding partial second derivatives. In our case, they can be obtained directly from equation (1). Given that each in-plane angle ${{\theta }_{ijk}}$ contributes energy change to two bonds (i-j and j-k) in the bond-order representation of graphene lattice, while each dihedral angle ${{\alpha }_{ijkl}}$ uniquely affects a single i-j bond, we get ${{k}_{bond}}=\frac{{{\partial }^{2}}{{V}_{ij}}}{\partial {{r}^{2}}}$, ${{k}_{ang}}=2\frac{{{\partial }^{2}}{{V}_{ij}}}{\partial \;{\rm cos} \;{{\theta }^{2}}}$, and ${{k}_{dih}}=\frac{{{\partial }^{2}}{{V}_{ij}}}{\partial {{\alpha }^{2}}}.$ All of the derivatives are evaluated at the corresponding prescribed values of the distances and angles. In the case of equation (1), these constants can be readily calculated analytically. From the mathematical standpoint, any quantitative agreement (in terms of the atomic positions as functions of time) between a system described by equation (1) and that described by equation (2) with the constraints defined above is at best expected at infinitesimally small atomic displacements, i.e. at very low simulated temperatures. Nevertheless, the extremely low computational cost of utilizing equation (2) makes it attractive in large-scale multicomponent atomistic or coarse-grained modeling. Furthermore, fine-tuning the stiffness parameters to better fit the experimental or ab initio calculated structural properties of membranes may be possible to some degree by analytically expressing these properties as functions of sought parameters (as presented in [22], for instance).

With the above constraint stiffness values, equation (2) is an incomplete harmonic constraint approximation of equation (1); hence, we refer to equation (2) as a bond-order-informed set of harmonic constraints. As noted in [15], equation (2) includes no energy cross-terms, i.e. the important interdependence of the pair and angular contributions explicitly present in equation (1) is removed from the decoupled energy expression in equation (2). It is then noteworthy that equation (2) represents a less accurate approximation of the original bond-order potential than that provided by the Green's function representation [23], which is an exact harmonic approximation of the original interactions and not a set of relative local constraints. As further shown, decoupled spring-like constraints expectedly represent a membrane model that is more resistant to thermally induced deformation and thus result in a stiffer graphene layer compared with the 'parent' bond-order potential.

Our goal here is not to develop a set of finely tuned harmonic constraints but a dynamic test for the interactions described by equation (2) with a set of constants determined directly from the bond-order representation. Therefore, we used the optimized second-generation Brenner formulation [6] to obtain ${{k}_{bond}}=698.129\;{\rm N}\;{{{\rm m}}^{-1}},$ ${{k}_{ang}}=8.080\;{\rm eV}$, and ${{k}_{dih}}=0.360\;{\rm eV}\;{\rm ra}{{{\rm d}}^{-2}}$, resulting in a 0 K bending rigidity of 2.2 eV, evaluated analytically, as described elsewhere [22]. Although it is somewhat higher than the value of 1.4 eV obtained from the second-generation bond-order potential [3], one must keep in mind that no consensus on the value of graphene's bending rigidity currently exists [24]. Moreover, the frequently cited experimental value of 1.2 eV was obtained indirectly from the phonon spectra of graphite [25]. A more refined set of constants can likely be obtained from any other suitable choice of original interactions, as suggested earlier. The validation procedure and the discussion below are generally applicable to any interatomic interaction potential that is properly approximated by equation (2).

Thermally induced dynamic rippling in atomically thin membranes is a naturally occurring phenomenon, arising from the phonons' fundamental inability to exist in two strict dimensions [26]. This thermally driven effect was theoretically quantified previously [16, 18] and was shown to generally correspond (for some nanoscale membrane dimensions, at least) to the harmonic theory of thermally fluctuating membranes [27]. We therefore selected the overall out-of-plane rippling magnitudes, the wave spectral distributions and the normal-normal correlations (as defined in [16] and [27]) as our quantities of interest. As an added point of comparison, we also simulated graphene with the use of a second-neighbor MEAM, as parameterized elsewhere [9].

The main simulated system is similar to that used in [18] for manipulating the dynamic rippling morphology in graphene (figure 1(a)). A suspended rectangular graphene membrane (with dimensions of 13.6 nm × 15.8 nm (8192 atoms)) was simulated dynamically for 1.0 ns with a time step of 1 fs. Prior to the production simulations, all of the graphene samples were statically pre-relaxed to a maximum interatomic force of 2.0 fN with use of the respective interatomic potential, resulting, in all cases, in a flat-phase graphene C-C bond length of 1.420 ${\rm {\ \mathring{\rm A} }}$. Because long vibrational modes tend to have a long survival time in MD simulations despite stringent thermal control, a decorrelated state in terms of the initial out of plane velocity components is important for simulating dynamic rippling in atomically thin membranes. Therefore, the initial Maxwell–Boltzmann-distributed atomic velocities were generated with the use of a high-quality random number generator [28]. All of the simulations were performed in an NVT ensemble in which a Nosé–Hoover thermostat [29] was utilized to maintain the prescribed temperatures. In order to minimize the direct effect of thermostatics on the dynamic ripples, thermal control was applied to the atoms within a 1 nm wide region along the membrane perimeter, as shown in figure 1(a), while the rest of the atoms were free of its direct effect. All of the quantities of interest were calculated for the atoms free of direct thermostatic control. All of the time averages were calculated within the last 0.75 ns of the simulated time, allowing the system to reach a steady state during the first 0.25 ns. Two smaller membranes were simulated: 11.1 nm × 12.8 nm (5408 atoms) and 5.5 nm × 6.4 nm (1352 atoms). In this work, the effective membrane dimensions were kept below the Ginzburg criterion limit of 20 nm, beyond which the anharmonic effects are expected to dominate the rippling process [16] and make the comparison with the theory questionable.

Zoom In Zoom Out Reset image size

According to the harmonic response theory of thermally fluctuating membranes, the out-of-plane mean square deviation (MSD) of thermal ripples in a thin membrane is $\left\langle {{h}^{2}} \right\rangle \propto TA/\kappa$, where T, A and $\kappa$ are its temperature, surface area and bending rigidity, respectively [16, 27]. It is important to note that the term 'harmonic' in the theory of fluctuating membranes refers to small fluctuations decoupled between the in-plane and the out-of-plane directions (see the supporting information for [16]). Such an approximation can be mathematically valid for infinitesimally small ripples in systems described by expressions such as equation (1). However, it is different from the decoupling presented by equation (2) compared with equation (1), in which the decoupling is between the local linear and angular energy terms both in-plane and out-of-plane. In the sense of the theory, both equations (1) and (2) possess the capacity for strongly anharmonic behavior at sufficiently high temperatures and large membrane dimensions.

Shown in figure 1(b) is a set of time-averaged out-of-plane MSD ${{\left. \left\langle {{h}^{2}} \right\rangle \right|}_{T}}$, curves for the interatomic interactions described by equation (2), the bond-order potential and the MEAM. Both the bond-order potential and the MEAM demonstrate expected behavior for the simulated membrane, although the MEAM appears to present a considerably stiffer model in terms of the ripple heights. The harmonic approximation yields a less than ideally linear dependence of 〈${{h}^{2}}$〉 on the temperature and is only in agreement with the 'parent' bond-order model at the lowest temperature of 50 K, as expected. At higher temperatures, the differences between equations (1) and (2) become quite considerable. The differences between the root mean square (RMS) displacements $\bar{h}=\sqrt{\left\langle {{h}^{2}} \right\rangle }$ obtained from the harmonic model and the bond-order model at the limiting temperatures of 50 K and 700 K are 13% and 200%, respectively. Importantly, at the low temperature of 1.0 K, the values of $\bar{h}$ are $4.01\times {{10}^{-2}}\;{\rm { \mathring{\rm A} }}$ and $3.91\times {{10}^{-2}}\;{\rm { \mathring{\rm A}}}$ for the bond-order and the harmonic model, respectively, and are in excellent agreement, as predicted above. From figure 1(b), the values of $\bar{h}\;$at 300 K are $0.44\;{\rm { \mathring{\rm A} }}$, $0.79\;{\rm { \mathring{\rm A} }}$, and $0.27\;{\rm { \mathring{\rm A} }}$, as obtained with the harmonic constraint model, the bond-order potential and the MEAM, respectively, while the C-C bond length in the flat phase of graphene is $1.42\;{\rm { \mathring{\rm A} }}$. As can be seen, the MEAM effectively represents the stiffest model of the three presented here. For general reference, the statistically characteristic ripple height in a sample of comparable size, as obtained from Monte Carlo simulations with the use of a LCBOPII bond-order potential in [16] at room temperature, is $0.7\;{\rm { \mathring{\rm A} }}$. As stated earlier, the bending rigidity of graphene (including its dependency on sample size and temperature), and the actual amount of thermal rippling and its spatial distribution are pending further experimental investigation.

To date, the only direct measurement of dynamic ripples in graphene was reported by P Xu et al [19], as obtained from the scanning tunneling microscopy (STM) measurements. The measured ripple root-mean-square deviations (RMSDs) reported in [19] are about an order of magnitude larger than those calculated previously [16–18] and those presented here. This level of discrepancy between the simulation and experiment can be expected due to scaling of the rippling magnitudes with the membrane size [17, 27]; thus, a direct quantitative comparison is only valid for membranes of equal size. Moreover, the STM tip-sample voltage can somewhat affect the measured data compared with an ideally free-standing membrane. It is worth mentioning that the experimentally measurable Debye–Waller factor that determines the effect of thermally induced lattice vibrations on the results of scattering experiments appears to have limited value in providing insight into the dynamic corrugation in the case of atomically thin layers, which are highly anisotropic in- and out-of-plane. Theoretical calculations performed within the harmonic approximation, which are used in determining the Debye–Waller factor, yield 0 K atomic MSDs of $1.6\;\times {{10}^{-3}}\;{{{\rm { \mathring{\rm A} }}}^{2}}$ and $4\times {{10}^{-3}}\;{{{\rm { \mathring{\rm A} }}}^{2}}$ in the in- and out-of-plane directions, respectively [30]. Later theoretical calculations yield a 0 K value of $2.2\times {{10}^{-3}}\;{{{\rm { \mathring{\rm A} }}}^{2}}$ of the overall MSD [31]. For comparison, at 1 K, both the bond-order and the harmonic constraints yield $\left\langle {{h}^{2}} \right\rangle =1.6\;\times {{10}^{-3}}\;{{{\rm { \mathring{\rm A} }}}^{2}}$, in accord with the results in [30]. The reasonable agreement near 0 K, however, should not suggest the use of the Debye–Waller factor as an accurate measure of the dynamic corrugation at higher temperatures. At 300 K, well below the Debye temperature, the corresponding atomic MSD estimated from the TEM studies of free-standing graphene is $4.4\times {{10}^{-3}}\;{{{\rm { \mathring{\rm A}}}}^{2}}$ [31], which is at least an order of magnitude below the amount of simulated out-of-plane rippling presented in figure 1(b) and in [16, 17]. Furthermore, this TEM result is about four orders of magnitude below the STM-measured $\left\langle {{h}^{2}}\; \right\rangle$ [19]. The cause for such a significant discrepancy is likely the lack of sensitivity of the TEM dark-field imaging technique on the relatively slow-moving long flexural waves responsible for the large out-of-plane ripples observed in the simulations and the STM measurements. The spectral distribution of the out-of-plane ripples therefore deserves further discussion.

The morphology of thermally fluctuating membranes can be further quantified with a spectral distribution of the normal-normal correlation. Here, the vectors that are normal to the surface are calculated locally from the atomic position snapshots taken throughout the steady-state simulation process at 300 K. As shown in figure 2, there is very good agreement between all three interatomic interactions, including the Bragg peak positions at ${{q}_{1}}=2.95\;{{{\rm { \mathring{\rm A} }}}^{-1}}$ and ${{q}_{2}}=5.9\;{{{\rm { \mathring{\rm A} }}}^{-1}}$, which is also consistent with [16]. For longer wavelengths ($0.1\;{{{\rm { \mathring{\rm A} }}}^{-1}}$ to $1\;{{{\rm { \mathring{\rm A} }}}^{-1}}$), qualitative trends in agreement with the theory of thermally fluctuating membranes [27] are also observed for all models in question. In addition to the statistics of the normal vectors, the spectral content of the ripples can be further quantified by calculating the Fourier spectral density of the ripple heights. Shown in figure 3 is a set of spectra obtained at 300 K for the models in question. The results are qualitatively similar, although the bond-order model reveals a considerably higher amount of long-wavelength rippling $\left( q\lt 0.1\;{{{\rm { \mathring{\rm A}}}}^{-1}} \right)$, which is apparently responsible for the overall higher amounts of rippling observed in figure 1(b) for the bond-order model. The particular peaks in figure 3 are shown for general comparison with the results in [16], which revealed a characteristic wavelength of ∼80 ${\rm { \mathring{\rm A} }}$ for a membrane of similar size.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In general, one expects the out-of-plane relaxation mechanisms to be different for various interatomic interactions. For instance, purely additive dihedral geometry constraints as part of implementing bending rigidity in equation (2) can have an interesting effect on relaxation in the dynamic case. From basic calculus, for any spectral component of amplitude ${{h}_{0}}(q)$, the mean square deviation of the local normal from the 001 direction is proportional to $\left( h_{0}^{2}(q){{q}^{2}} \right)\propto \;\left( h_{0}^{2}(q)/{{\lambda }^{2}} \right)$. This average can also be seen as a measure of the membrane's dihedral-based instantaneous roughness for a given Fourier component. Reducing the overall statistical average of the dihedral angles (as a way to satisfy the dihedral constraint) can therefore be achieved by suppressing the rippling amplitudes as well as by introducing longer waves of possibly larger amplitudes, as long as the overall $h_{0}^{2}/{{\lambda }^{2}}$ is reduced, even if $h_{0}^{2}\;$ is increased [32]. The relative contributions from these two relaxation mechanisms are determined by the amount of anharmonic effects present in the interatomic interactions [16]. Let us consider figure 4, in which we show the relevant portions of the ripple spectra for the harmonic model with and without the dihedral constraint. An increasing presence of the long flexural modes is observed with the dihedral constraint, as shown. The longer flexural waves clearly require compatibility with the membrane dimensions, representing a redistribution of deformation energy in an isothermal system. Given the discussion above, an overall increase in $\left\langle {{h}^{2}} \right\rangle$ does not imply a more flexible membrane governed by equation (2). In fact, it implies the opposite, as shown in figure 5, in which we plot the RMS value of the out-of-plane angle and the out-of-plane RMSD $\bar{h}$ as functions of varying ${{k}_{dih}}$. An increase in the dihedral constraint stiffness can thus result in increasing $\bar{h}\;$ and decreasing the RMS values of the dihedral angles. This effect must be taken into account when fine-tuning the harmonic constraint stiffness set. We note that the bond-order model (which corresponds to ${{k}_{dih}}=0.360\;\;{\rm eV}\;{\rm ra}{{{\rm d}}^{-2}}$) yields an RMS dihedral deviation of 0.138 rad, which is compared with 0.163 rad in figure 5, while the MEAM yields 0.111 rad, resulting once again in the stiffest model both in terms of the dihedral RMS and the rippling magnitudes (see figure 1(b)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In general, we expect the differences between the harmonic constraint model and the original bond order representation to somewhat decrease with the decreasing membrane size. All of the out-of-plane displacement components that affect the bond-angle coupling in equation (1) are expected to decrease [27], somewhat reducing the overall differences between the representations provided by equations (1) and (2). See figure 6 for a comparison of the out of-plane MSD as a function of temperature for smaller membranes. Although the differences between the harmonic model and the bond-order potential are still considerable, we note that the fine-tuning of the constraint stiffness values, along with the use of thermal over-relaxation (e.g. with use of a Langevin thermostat [33]) if appropriate, may bring the two models even closer for membranes of various dimensions. Therefore, given the wide availability of the harmonic constraint model in various simulation packages, along with its computational efficiency, the presented parameterization may be an attractive trade-off between the accuracy and simplicity of the representation.

Zoom In Zoom Out Reset image size

We presented a complete set of harmonic bond, angle and dihedral angle constraints, as calculated directly from the bond-order potential, for a simplified model of free-standing graphene that is suitable for a variety of existing molecular mechanics simulation packages. We performed a number of carefully designed atomistic molecular dynamics simulations and quantified the dynamic thermal rippling process in free-standing graphene. The described studies present a way to assess the dynamic performance of the harmonic constraint model at various finite temperatures, in addition to the structural validation presented in [15]. Our results indicate that a simple harmonic constraint model with full decoupling between the energy terms associated with linear and angular lattice perturbations yields a stiffer graphene membrane in comparison with the bond-order potential it is based on. Nevertheless, the rippling averages are generally consistent with previous reports. In addition, the normal-normal correlations and height distributions of ripples obtained with the harmonic constraint model are in good agreement with the bond-order model and the MEAM. Importantly, we demonstrate that increasing the dihedral angle constraint results in dynamic membrane relaxation via long flexural waves of increased amplitude. This result appears to be counterintuitive but is in agreement with the theory of fluctuating membranes and thus should be taken into account when parameterizing a set of harmonic constraints for atomically thin layers.

The harmonic constraint model is widely available in a number of existing biomolecular simulation software packages. Therefore, given its overall reasonable dynamic accuracy and very low computational cost, we believe that such a model may be suitable for use in multicomponent dynamic simulations, especially in systems that include atomically thin membranes in an aqueous environment.

The author thanks A N Chiaramonti for the enlightening discussion on the TEM dark-field imaging technique. I am also grateful to A F Kazakov and V K Tewary for their numerous useful discussions. This work is a contribution of the National Institute of Standards and Technology, an agency of the US government; this work is not subject to copyright in the USA.