Surface energy is a fundamental material property that determines important functions such as catalytic, sensing, and imaging properties. Over the past century, various experimental studies and models including the broken bond theory and Wulff construction have been developed to analyze surface free energies. However, it remains a challenge to measure or predict thermal fluctuation effects on surface energies. In particular, crystals of functionalized building blocks, such as self-assembling proteins and DNA-functionalized nanoparticles, assembled via the specific surface interactions of the building blocks, are highly sensitive to thermal fluctuations. In the case of DNA-functionalized nanoparticles, it has been shown that the crystals are formed as a result of thermally active hybridizations. We show here that the surface energy along different planes can be obtained from the ratio of hybridization events. The surface energy fluctuations in these systems are shown to bear a nearly linear correlation with the fluctuations in DNA hybridization events in the bulk. We further demonstrate that short DNA chains and high DNA loading increase the volume density of the DNA sticky ends. The relationship between thermally active hybridizations and surface energy found here can be used to aid the design of single crystals of functionalized colloids with active surface groups.

## INTRODUCTION

surface energy has been a topic of extensive study since the seminal work by Gibbs.^{1} As one of the most fundamental properties of a solid, surface energy determines the equilibrium shape of a crystal and plays an essential role in faceting and roughening, as well as in crystal growth phenomena.^{2} The fluctuations of surface energy, while important for growing a stable and sharp single crystal,^{2,3} have seldom been reported or studied. This is because the accuracy and precision of most experimental methods fall short producing relative errors ranging from 5% to 30%.^{4} Unfortunately, this difficulty only increases when trying to measure thermal fluctuations at a surface. On the theory side, density functional theory has helped to establish a database of surface energies,^{5} but energy fluctuations have not yet been reported using such an approach.

With the success of the first single crystal growth composed of DNA-functionalized nanoparticles (DNA-NP), these man-made “programmable atom equivalents” provide a robust system for growing faceted single crystals in soft materials.^{6} Over the past two decades, this elegant strategy of using complementary DNA linkers to guide the assembly of nanoparticles has increased our ability to control crystal structure, resulting in a transition from disorganized amorphous products^{7,8} to organized crystalline symmetries,^{9–13} as well as micron-scale single crystals with fewer defects.^{6} Excitingly, the DNA-NP system provides us with a unique opportunity to explore the origin of surface energy fluctuations, by probing the thermodynamics relevant to DNA chains via molecular dynamics (MD) simulations. Using DNA hybridizations to mimic atomic-scale chemical bonds, one can track the details of a single DNA attachment and detachment event.^{14–18} Using MD simulations, one is able to capture properties associated with thermal fluctuations, providing a deeper understanding of the atomic crystal systems as well.

It is widely accepted that surfaces bearing large thermal fluctuations do not benefit the stability of single crystal polyhedrons.^{2,19,20} Therefore, where do these surface energy fluctuations come from and how can we design a system that minimizes their effect whilst optimizing crystal growth? In this work, we apply a molecular model with robust MD simulations to extract a mathematical relationship bridging design parameters (e.g., length of an individual DNA chain, DNA coverage on each NP, NP size) and surface energy fluctuations. Based on this relationship, we discuss optimal design parameters that minimize these fluctuations for the design of single crystal polyhedra.

## MODEL

Herein we used a scale-accurate coarse-grained model with explicit DNA chains, which faithfully captures relevant sizes and stiffnesses for the double-stranded (dsDNA) and the single-stranded DNA parts in the design. The attraction between complementary DNA linkers (or sticky ends) is modeled by the Lennard-Jones potential. This model excellently represents DNA hybridization kinetics and the energy level per DNA hybridization, and it has been reported to be robust and reliable in accurately recovering the self-assembly processes of various symmetries with DNA-NPs.^{15,16}

The quantity of interest is the surface excess free energy per unit area of a specific crystal facet *γ*. By custom, *γ* is referred to as the “surface energy.” A direct measurement of *γ* requires to separate a bulk crystal into two parts along a particular plane. In the modeled system, periodic boundary condition along the z-axis was removed from the bulk crystal to expose the facet of interest on two sides (Figure 1). To calculate *γ* of the exposed facet, we measured the energy difference between the bulk system and the system with two exposed surfaces, then divided the value by twice the surface area. To keep within the capability of molecular modeling, several assumptions were made regarding entropic contributions. In particular, a perfect crystalline structure and constant lattice parameter are assumed.^{6} That is, we are neglecting the entropic contributions from the lattice vibrations. *γ* is computed taking the time average of the potential energy.

The Wulff construction is an elegant method to determine the equilibrium crystal shape of fixed volume, known as the Wulff shape,^{23,24} if the orientational dependence of the surface free energy is known. For each plane $ h k l $, plotting a point in the $ h k l $ direction at a distance which is equal to the associated surface energy *γ*(*hkl*) gives the Wulff plot. Planes perpendicular to the radial direction are drawn on each point of the Wulff plot. The inner envelope of all these planes gives the equilibrium shape, or the Wulff shape. With surface energy values estimated from MD simulations, we apply here the WulffMaker software^{25} to calculate and display the corresponding Wulff polyhedron.

## RESULTS

Wulff construction indicates that, instead of the absolute surface energy values, it is the surface energy ratios of different facets that determines the equilibrium shape of a crystal. For example, the body-centered cubic (BCC) metal should exhibit a ratio of *γ*(110):*γ*(100):*γ*(111) = 1:1.41:1.22^{26,27} corresponding to a rhombic dodecahedron; the face-centered cubic (FCC) metal should exhibit *γ*(111):*γ*(100):*γ*(110) = 1:1.15:1.22, which corresponds to a truncated octahedron. Naturally, surface energy fluctuations propagate to the *error bars* of this ratio. In this section, we derive a correlation between the fluctuations on surface energy ratios and fluctuations of DNA hybridization events with several assumptions made. We performed simulations on distinct systems to verify the model assumptions and to validate the obtained mathematical relationships. We then bridge the fluctuations of DNA hybridizations with design parameters directly controllable in experiments. As such, we found the relationship between design parameters and surface energy ratios.

### surface energy fluctuations

Due to the expandable nature of spherical nucleic acid,^{28} building blocks only interact with their first nearest neighbors^{12} thus low-index planes such as (100), (110), and (111) become cusp orientations in the Wulff plot,^{26} determining the equilibrium shape. In general, the surface energy of any facet (*hkl*) in a DNA-NP superlattice is calculated by Eq. (1),

and

where *E _{hkl}* and

*E*

_{bulk}are the energies of a system with two (

*hkl*) surfaces exposed and the energy of the bulk system, respectively.

*p*is the percentage of hybridized DNA pairs,

_{H}*N*

_{DNA pairs}is the maximum number of DNA hybridizations that could be formed, and

*E*

_{hybridization}is the hybridization energy per linker attachment (e.g.,

*E*

_{hybridization}= − 42.3

*k*J/mol for DNA linker sequence −TTCCTT

^{16}). $ p H $ is the ensemble average of

*p*and Δ

_{H}*p*is the standard deviation of

_{H}*p*defined by $\Delta p H = p H 2 \u2212 p H 2 1 / 2 $.

_{H}By propagation of uncertainty (e.g., for $f=c ( A + B ) ,\Delta f=c \Delta A 2 + \Delta B 2 $, if *A* and *B* are uncorrelated.), $\Delta \gamma ( h k l ) = \eta h k l \xd7 \Delta p H 2 ( h k l ) + \Delta p H 2 ( bulk ) $ as the two measurements of *p _{H}* are independent from each other.

As the bulk system and the (*hkl*) system have a similar system size, we can approximate that $\Delta p H 2 ( h k l ) \u2248\Delta p H 2 ( bulk ) \u2248\Delta p H 2 $. We further define $ k h k l = p H ( h k l ) p H ( bulk ) $, which is a mean field value consistent with the classical broken bond model.^{26,27} It is thus reasonable to consider *k _{hkl}* to be identical for spherical nanoparticles with different design parameters, as long as the crystalline systems bear the same dimension (i.e., the number of unit cells on each dimension of the system box keeps the same). As such, the relative surface energy fluctuation of facet (

*hkl*) can be approximated as

Thus, surface energy fluctuations originate from the fluctuations of DNA hybridization events. This is consistent with the thermally active hybridization crystallization process, which has been shown to be crucial for the successful crystallization of DNA-NPs.^{16}

### Relative error in surface energy ratio

surface energy ratio is given respective to the reference facet (ref.) which has minimum surface energy of all facets. For example, facet (110), the closest-packed plane in BCC, is used as the reference for BCC symmetry; similarly, facet (111) is the reference for FCC symmetry. surface energy ratio regarding any facet (*hkl*) is expressed as $ \gamma ( h k l ) \gamma ( ref. ) $, here after refer to as *g ^{hkl}*.

The uncertainty in *g ^{hkl}* is estimated by

and

(Propagation of uncertainty: for $f= A B ,\u2009 \Delta f f = \Delta A 2 A 2 + \Delta B 2 B 2 $, if *A* and *B* are uncorrelated.)

With the above derivation, the relative error in *g ^{hkl}* is seen to bear an approximately linear correlation with the relative error in

*p*(bulk) with coefficient

_{H}*C*.

^{hkl}For the sake of comparability and simplicity, all of the systems tested below were constructed into the BCC symmetry. Nevertheless, the protocol we illustrate and the phenomena we discuss in this work is applicable to any crystalline symmetry in the DNA-NP system.

As shown in Figure 2, different BCC lattice were constructed to study different facets of interest, and all are composed of 96 building blocks. Since we have noted different values of $ k h k l = p H ( h k l ) p H ( bulk ) $ for different system geometries, we provide an estimation for *k*_{100}, *k*_{110}, and *k*_{111}, as well as the coefficients *C*^{100} and *C*^{111} for the specific geometries study here. To better evaluate the model, these theoretically estimated parameters are validated in the section titled Validation by simulations by straightforward simulations.

Consistent with the broken bond model,^{26,27} a number of building block connections are cut off when the bulk material is cleaved to expose a specific facet. As each building block’s nearest neighbors consist of eight building blocks with the complementary DNA linker type in the BCC lattice, we approximated that each neighbor contributes one eighth of its total connections. Note that building blocks do not connect with their second nearest neighbors in this DNA-NP design. We counted the number of neighbors cut off to estimate the coefficients as in Eq. (4). Parameters of 1/4, 1/2, and 3/8 are simply the portions of neighbors got disconnected by exposing the specific surface plane,

Using Eq. (3), the coefficient for the BCC (100) system is calculated to be *C*^{100} = 14.14, and for (111) system it is *C*^{111} =16.00.

### Validation by simulations

Until now the mathematical relationship (Eq. (3)) has been derived under several assumptions. In order to verify this formula, we simulate five distinct DNA-NP designs with different design parameters and directly calculate their surface energy values. For each DNA-NP design, we carry out MD simulation runs on approximately five independently generated initial configurations with identical design parameters (DNA chains are randomly attached on each NP). Simulation statistics are summarized in Table I. The error bars in the table are calculated using the mean of the standard deviations of the five identical simulation runs.

. | Design parameters^{a}
. | . | . | . | ||
---|---|---|---|---|---|---|

No. . | NP size (nm) . | DNA coat . | dsDNA (bp) . | ρ^{b} (nm^{−3})
. | p_{H}^{c} (%)
. | Δp/_{H}p_{H}^{c}
. |

1 | 10 | 60 | 18 | 0.017 | 37.84 | 0.0156 |

2 | 10 | 60 | 58 | 0.0014 | 13.23 | 0.0368 |

3 | 10 | 90 | 18 | 0.025 | 41.56 | 0.0110 |

4 | 10 | 90 | 58 | 0.0021 | 16.84 | 0.0254 |

5 | 20 | 150 | 18 | 0.023 | 32.90 | 0.0105 |

. | Design parameters^{a}
. | . | . | . | ||
---|---|---|---|---|---|---|

No. . | NP size (nm) . | DNA coat . | dsDNA (bp) . | ρ^{b} (nm^{−3})
. | p_{H}^{c} (%)
. | Δp/_{H}p_{H}^{c}
. |

1 | 10 | 60 | 18 | 0.017 | 37.84 | 0.0156 |

2 | 10 | 60 | 58 | 0.0014 | 13.23 | 0.0368 |

3 | 10 | 90 | 18 | 0.025 | 41.56 | 0.0110 |

4 | 10 | 90 | 58 | 0.0021 | 16.84 | 0.0254 |

5 | 20 | 150 | 18 | 0.023 | 32.90 | 0.0105 |

. | γ(100)^{d}
. | γ(110)^{d}
. | γ(111)^{d}
. | Ratios of surface energies^{d}
. |
---|---|---|---|---|

No. . | (mJ m^{−2})
. | (mJ m^{−2})
. | (mJ m^{−2})
. | γ(110) : γ(100) : γ(111)
. |

1 | 0.531 ± 0.077 | 0.354 ± 0.073 | 0.440 ± 0.089 | 1:1.50 ± 0.38:1.24 ± 0.36 |

2 | 0.057 ± 0.018 | 0.039 ± 0.017 | 0.049 ± 0.021 | 1:1.45 ± 0.79:1.25 ± 0.77 |

3 | 0.853 ± 0.082 | 0.573 ± 0.078 | 0.716 ± 0.095 | 1:1.49 ± 0.25:1.25 ± 0.24 |

4 | 0.103 ± 0.022 | 0.070 ± 0.021 | 0.087 ± 0.026 | 1:1.47 ± 0.55:1.24 ± 0.53 |

5 | 0.628 ± 0.058 | 0.434 ± 0.055 | 0.531 ± 0.068 | 1:1.45 ± 0.23:1.22 ± 0.22 |

. | γ(100)^{d}
. | γ(110)^{d}
. | γ(111)^{d}
. | Ratios of surface energies^{d}
. |
---|---|---|---|---|

No. . | (mJ m^{−2})
. | (mJ m^{−2})
. | (mJ m^{−2})
. | γ(110) : γ(100) : γ(111)
. |

1 | 0.531 ± 0.077 | 0.354 ± 0.073 | 0.440 ± 0.089 | 1:1.50 ± 0.38:1.24 ± 0.36 |

2 | 0.057 ± 0.018 | 0.039 ± 0.017 | 0.049 ± 0.021 | 1:1.45 ± 0.79:1.25 ± 0.77 |

3 | 0.853 ± 0.082 | 0.573 ± 0.078 | 0.716 ± 0.095 | 1:1.49 ± 0.25:1.25 ± 0.24 |

4 | 0.103 ± 0.022 | 0.070 ± 0.021 | 0.087 ± 0.026 | 1:1.47 ± 0.55:1.24 ± 0.53 |

5 | 0.628 ± 0.058 | 0.434 ± 0.055 | 0.531 ± 0.068 | 1:1.45 ± 0.23:1.22 ± 0.22 |

^{a}

For example, each building block consists of a nanoparticle sized at 10 nm, coated by 60 DNA chains, and each DNA chain contains 18 base pairs (bp) for the dsDNA part.

^{b}

Volume density of DNA linkers, introduced in Eq. (5).

^{c}

*p _{H}* and Δ

*p*/

_{H}*p*refer to the value in bulk system, calculated from the DNA hybridization potential energy.

_{H}^{d}

Surface energy values correspond to the 96 building block system discussed in the text. According to the central limit theorem, standard deviation decreases with sample size as $\sigma \u221d 1 N $, where *σ* is the standard deviation and *N* is the number of building blocks in the simulation system.

Note that the surface energy values shown in Table I correspond to the 96 building block systems in Figure 2. According to the central limit theorem, standard deviation decreases with the number of building blocks: $ \Delta g h k l g h k l \u221d 1 N $, where *N* is the number of building blocks in the simulation system. Although *N* = 96 was the largest system we could simulate due to memory budget, the critical size for a single crystal to retain a stable structure is *N* ∼ 10 000.^{19} Structural changes were reported for clusters below a few thousand atoms. Importantly, we should note that in the range of *N* ≈ 10 000, the fluctuation levels of the surface energy ratios are about one tenth of the values listed in the table.

Regression of $ \Delta g h k l g h k l $ on $ \Delta p H p H $ fit excellently into a linear model as shown in Figure 3, and the coefficients highly agree with the derived values. This statistical analysis has verified the linear correlation between the relative errors in the surface energy ratios and the relative fluctuation of DNA hybridizations.

We have confirmed that the minimization of DNA hybridization fluctuations is essential to minimize the errors in surface energy ratios. We next predict DNA hybridization fluctuations on design parameters controllable in experiments. Accordingly, bulk systems with numerous sets of design parameters are simulated, including NP sizes at 5, 10, 15 nm, DNA chains composed of 18, 38, 58, 78 base pairs and various levels of DNA coats.

Based on tests with several variables, we find that the volume density of DNA linkers *ρ* (Eq. (5)) is an excellent indicator to predict Δ*p _{H}*/

*p*,

_{H} where *V*_{layer} = *V*_{max} − *V*_{bb}, $ V max = 4 3 \pi R max 3 $, and $ V bb = 4 3 \pi R 3 $.

In Eq. (5), *R* represents the hydrodynamic radius of a DNA-NP building block,^{12} *R*_{max} represents the maximum radius that a fully stretched building block can reach^{12} and *r* is the number of DNA chains coated on each nanoparticle. Correspondingly, *V*_{bb} is the hydrodynamic volume of a DNA-NP building block, *V*_{max} is the maximum volume a building block can reach with DNA chains fully stretched out, and *V*_{layer} is the estimated volume where sticky ends appear.

As shown in Figure 4, *ρ* and Δ*p _{H}*/

*p*fits into a power law model $log ( r \xd7 \Delta p H / p H ) =\u22120.35\xd7log ( \rho ) \u22121.54$, with the coefficient of determination

_{H}*R*

^{2}= 0.983. Taking the cue from this analysis, high volume density of sticky ends is preferred for small errors in surface energy ratios. From Eq. (5), high DNA coat, short DNA chains as well as large NP if DNA loading density retain constant will benefit a high volume density of sticky ends. Interestingly, experiments fail to grow single crystal Wulff shapes with 20-nm NPs coated by long DNA chains (e.g., 58 or 98 dsDNA pairs).

^{6}

We applied the WulffMaker^{25} to display Wulff polyhedra for surface energy ratios with different fluctuation levels, as shown in Figure 5. Ideally, the Wulff shape for a BCC metal is the rhombic dodecahedron. We note that in Figure 5, we are not constructing Wulff shapes with fluctuations in surface energy since Wulff shapes are constructed from ratios of mean surface energies. In the hypothetical case that the surface energies lower bounds could be accessible experimentally by growing Wulff shapes while suppressing the fluctuations (say with different solvent conditions and/or grown at low temperatures), and then the Wulff shapes are allowed to relax say by annealing them, the relaxed Wulff shape will become more globular.

## SIMULATION DETAILS

The MD simulations were performed with the HOOMD-Blue package^{29,30} (available at http://codeblue.umich.edu/hoomd-blue) on graphics cards. Constant volume and temperature ensemble (NVT) was applied with the temperature controlled via a Nosé-Hoover thermostat.^{31,32} All the simulations were performed at the same temperature, providing a universal reference state. Additional details of the model and the simulation protocol can be found elsewhere.^{6,15,16} The molecular model and the simulation protocol have been shown to capture the observed crystal symmetries (BCC, FCC, CsCl, AlB_{2}, Cr_{3}Si, Cs_{6}C_{60}).^{15,16} In the present study, since a crystal composed of ∼100 nanoparticles may not stay stable, the positions of the nanoparticle centers are fixed while performing MD simulations yet the particles are free to rotate.^{6} To keep the system within the memory limitation of a Tesla K40 graphics card, 96 independently generated DNA-NPs were simulated for each design with DNA chains randomly attached on each NP. Nonetheless, the system is larger than the critical nuclei size which is reported to be ∼40 NPs for the DNA-NP system.^{33}

## CONCLUDING REMARKS

This work provides a method for understanding and predicting surface energy fluctuations in DNA-NP systems that originate from thermally active hybridization events, while the latter is found to be easily predicted using the *volume density of DNA sticky ends*. Indeed, these simulation results are the first quantitative measurements that demonstrate the relationship between thermal fluctuations at surfaces and the thermal fluctuations of DNA hybridization events. We found that shorter DNA chains and higher loading of DNA chains on NPs are the main factors to reduce the errors in surface energy ratios. Our studies show that the larger the number of particles in a single crystal, the less surface energy fluctuation errors. Therefore, finding ways to grow larger crystals is essential to develop single crystals of nanoparticles with specific optical and electronic properties. Furthermore, we provide a general method to understand surface energy fluctuations. This method can be adapted to describe the directed-assembly of anisotropic particle shapes^{11,34} and colloids covered by multiple types of ligands.^{35} It would be interesting to determine how these fluctuations could be introduced in other level of coarse-grained models such as those in Refs. 36 and 37.

## Acknowledgments

M.O.d.l.C. and T.I.N.G.L. acknowledge support from the Air Force Office of Scientific Research (AFOSR) Award Nos. FA9550-11-1-0275 and FA9550-10-1-0167. T.I.N.G.L. gratefully acknowledges the Ryan Fellowship and the Northwestern University International Institute for Nanotechnology.