We use molecular simulation to construct equilibrium phase diagrams for two recently introduced model materials with isotropic, soft-repulsive pair interactions designed to favor diamond and simple cubic lattice ground states, respectively, over a wide range of densities [Jain et al, Soft Matter9, 3866 (2013)]. We employ free energy based Monte Carlo simulation techniques to precisely trace the inter-crystal and fluid-crystal coexistence curves. We find that both model materials display rich polymorphic phase behavior featuring stable crystals corresponding to the target ground-state structures, as well as a variety of other crystalline (e.g., hexagonal and body-centered cubic) phases and multiple reentrant melting transitions.

Crystalline materials with simple cubic or diamond symmetries are often desired for photonic or other optical applications.1,2 Unfortunately, large-scale and inexpensive fabrication methods for such materials which allow one to prescribe both the physical dimensions and the symmetry of the periodically replicated features have been challenging to develop.3,4 Holographic lithography, where a photoresist film is exposed to an optical profile with the desired symmetry created by interference of multiple coherent laser beams, represents one promising top-down approach.5 On the other hand, self-assembly of micron or nanoscale particles from disordered fluid states into a targeted lattice structure is a widely studied bottom-up strategy.

Although practical challenges remain, computer simulations and experiments have shown that self-assembly of the low-coordinated lattice structures of interest is possible by selecting particles with specific anisotropic shapes (e.g., tetrapods, pyramids, etc.)6–9 or orientationally specific “patchy” interactions.10–17 What is perhaps surprising is the prediction from computer simulations that simple cubic and diamond crystalline phases can also result from self-assembly of particles with isotropic interactions18–30 including some that are purely repulsive.21–29 These latter systems are of practical interest because, while large quantities of nanoparticles with precise anisotropic interactions can be challenging and expensive to manufacture, those with approximately isotropic interactions (e.g, ligand-coated nanocrystals) can be readily synthesized. Furthermore, since such isotropic pair interactions are “effective” in nature, they can be tuned by, e.g., the choice of the ligand and the properties of the solvent.31–33 Finally, self-assembly of particles with isotropic interactions may be advantageous from a kinetic perspective when compared to patchy particles, since the energy landscapes of the latter can be particularly rugged.34 

Motivated by an insightful study on statistical mechanical design of isotropic interactions for low-coordinated ground states in two dimensions,35 we recently introduced an inverse approach28 for discovering simple pair potentials36 that exhibit three-dimensional diamond37 or simple cubic38 ground states stable over a wide range of densities. The qualitative characteristics of the interaction potentials we optimized (short-range, isotropic, convex, and repulsive) were inspired by the soft, entropic effective repulsions encountered in many polymer(ligand)-decorated colloidal(nano) particle systems. The optimized interactions we obtained for target ground states with diamond (φdia) and simple cubic (φsc) structures are only subtly different from one another. For a comparison, see insets (c) and (d) of Figure 1. In short, the former has a softer core and a slightly harder shoulder than the latter. As we discuss below, these differences have important implications for the crystalline phases favored at intermediate densities. We refer readers to our previous publication28 for details on the stochastic optimization routine used to obtain these potentials and the corresponding ground-state phase diagrams. A similar pair potential with a diamond ground state was independently discovered by Marcotte, Stillinger, and Torquato using an alternative inverse design strategy with different constraints,29 suggesting that the qualitative form of the optimized interaction is robust.

FIG. 1.

Pressure-temperature projection of the phase diagrams for the (a) diamond forming potential (φdia) and (b) the simple cubic forming potential (φsc). Points are from the reference phase coexistence simulations discussed in the text, and curves are phase coexistence traces obtained directly from the NPT-TEE MC simulations. Statistical uncertainties are smaller than the symbol size in all cases. Crystal phases are those described in the text and in Table I. Inset (c) presents both potentials36–38 and inset (d) elucidates the difference (φsc − φdia) as a function of interparticle separation r/σ.

FIG. 1.

Pressure-temperature projection of the phase diagrams for the (a) diamond forming potential (φdia) and (b) the simple cubic forming potential (φsc). Points are from the reference phase coexistence simulations discussed in the text, and curves are phase coexistence traces obtained directly from the NPT-TEE MC simulations. Statistical uncertainties are smaller than the symbol size in all cases. Crystal phases are those described in the text and in Table I. Inset (c) presents both potentials36–38 and inset (d) elucidates the difference (φsc − φdia) as a function of interparticle separation r/σ.

Close modal

In this Communication, we use free energy based Monte Carlo (MC) simulation methods to elucidate the effect of temperature (entropy) on the aforementioned diamond forming and simple cubic forming potentials28 by constructing their equilibrium phase diagrams. Despite the apparent simplicity of the interparticle interactions, we find that the models display rich, polymorphic phase behavior. Crystalline phases corresponding to the target simple cubic and diamond ground-state structures feature prominently on their respective phase diagrams as well as a variety of other interesting crystalline phases and multiple reentrant melting transitions.

In order to trace the loci of phase coexistence points in the temperature–pressure–density (TP − ρ) space, we first evaluate at least one reference inter-crystal and crystal-fluid coexistence point as follows. We locate crystal melting points along a specific isobar by computing the temperature at which the Gibbs free energy of the crystal equals that of the fluid. We employ Frenkel-Ladd39 and isothermal-isobaric temperature-expanded ensemble MC simulations40,41 to determine the temperature dependence of the Gibbs free energy of the crystalline phases and a combination of grand canonical transition matrix MC42 and isothermal-isobaric temperature expanded ensemble (NPT-TEE) MC40,41 simulations to determine the Gibbs free energy of the fluid phase as a function of temperature. The inter-crystal coexistence points are evaluated for a single low-T isotherm (T = 0.01ε/kB, where kB is the Boltzmann constant)—and other isotherms as needed—using a combination of Frenkel-Ladd39 and density expanded ensemble MC simulations.41 

We trace inter-crystal coexistence curves and fluid-crystal saturation curves by employing a recently introduced technique41 that also uses isothermal-isobaric temperature expanded ensemble (NPT-TEE) MC simulations, wherein the sub-ensembles are differentiated by a shift in temperature and pressure. In this method, one first completes a NPT-TEE simulation with a guess of the temperature-pressure relationship along the saturation line and later corrects this guess by employing histogram reweighting with the relevant probability distributions (e.g., volume and enthalpy) to identify the coexistence point associated with each sub-ensemble sampled. This approach allows us to compute coexistence points over a wide range of temperatures within a single simulation. For tracing inter-crystal coexistence curves, we use NPT-TEE simulations where temperature is the dominant variable and work with the resulting density probability distribution. We then use histogram reweighting to determine the pressure at which the chemical potentials of the two crystal phases are equal at each sub-ensemble temperature. For solving the fluid-crystal coexistence curves, we use another variant of NPT-TEE where pressure is the dominant variable. We utilize the enthalpy probability distribution in this simulation and compute the coexistence temperature at each sub-ensemble pressure by using histogram reweighting. This strategy requires a reference saturation state point and the corresponding Gibbs free energy, which we provide from the single isobar MC simulations explained above. We refer the reader to the earlier publication41 for a detailed description of this approach. The NPT-TEE simulation of the crystal phase is able to sample inter-crystal transitions or crystal-fluid transitions for cases when the specified pressure-temperature sub-ensemble relationship extends to conditions beyond the crystal's stability range. Thus, it often becomes apparent from results of the simulation if another crystal (which might not have been considered in the original free energy calculations) is more stable at a given state point—information that is then used to suggest new free energy calculations to refine the phase diagram. We also completed several quench simulations at different densities, wherein we allow a fluid (initially equilibrated at a high temperature of T = ε/kB) to relax at lower temperatures, and we verified that it assembles into the expected crystal (for more information, see the supplementary material43).

The structures we consider in the free-energy calculations consist of all lattices which occur as ground states for the optimized potentials,28 as well as other lattices which we found to be slightly less stable than the ground states (by an enthalpy difference of order 10−4ε). This composite list includes face-centered cubic (FCC), body-centered cubic (BCC), simple cubic (SC), diamond (DIA), hexagonal (H), A7, A20, βSn, and cI16 crystalline phases (Table I). We also allow for changes in box shape and in the lattice parameters during the simulations. The MC simulations typically contained 700–2000 particles, the actual number depending on the lattice-type. We simulated larger sized simulations with ⩾ 3000 particles and did not observe significant differences in the free energies.

Table I.

Lattice parameters for the equilibrium crystal phases of the optimized potentials, φdia and φsc (nomenclature is that of an earlier reference44).

Diamond forming potential model, φdia
Hexagonala c/a: [ 0.589 − 0.734] 
A7 b/a: 3.98, u: 0.317 
Hexagonal c/a: 1.4 
A20 b/a: 1.73, c/a: 0.622, y: 0.33 
Simple cubic forming potential model, φsc 
Hexagonala c/a: [ 0.845 − 0.89] 
A20 b/a: 2.27, c/a: 1.8, y: 0.9025 
βSn c/a: [ 0.66 − 0.7] 
Hexagonal c/a: [ 0.95 − 1.0] 
A20 b/a: 3.1, c/a: 1.02, y: 0.867 
Diamond forming potential model, φdia
Hexagonala c/a: [ 0.589 − 0.734] 
A7 b/a: 3.98, u: 0.317 
Hexagonal c/a: 1.4 
A20 b/a: 1.73, c/a: 0.622, y: 0.33 
Simple cubic forming potential model, φsc 
Hexagonala c/a: [ 0.845 − 0.89] 
A20 b/a: 2.27, c/a: 1.8, y: 0.9025 
βSn c/a: [ 0.66 − 0.7] 
Hexagonal c/a: [ 0.95 − 1.0] 
A20 b/a: 3.1, c/a: 1.02, y: 0.867 
a

In the lower density region.

We present the computed projections of the fluid-crystal and inter-crystal phase coexistence curves in the PT and T − ρ planes in Figures 1 and 2, respectively. First, we observe that the target phases (diamond for φdia and simple cubic for φsc) feature prominently on the corresponding phase diagrams. This is notable because such phases rarely occur in systems with isotropic interactions, and—when they have appeared in model systems—they generally display a narrow range of thermodynamic stability.28 In this case, however, the potentials were explicitly designed to exhibit the target structures in the ground state (at T = 0) over a wide range of densities.28 What our new calculations show is that the crystalline phases corresponding to the ground-state optimization targets also exhibit good thermal stability (i.e., comparable to other stable lattices on the phase diagram). This is encouraging because it suggests that focusing exclusively on the ground-state structure during optimization of the pair potential can be an effective strategy, even if the ultimate goal is to find interactions that trigger assembly into the target phase at finite temperature. We further note that several structures for which the potentials were not optimized (e.g., βSn, rhombohedral, and body-centered orthorhombic lattices), but that nonetheless appear in the ground-state phase diagram,28 display poor thermal stability with only βSn surviving to the lowest temperatures explored here. The ground-state optimization of the potentials for the target lattices against a large pool of competitive structures across a wide density range may well have had the secondary effect of increasing the relative thermal stability of the target phases. We may gain more insight into this hypothesis by constructing the equilibrium phase diagrams of interaction potentials designed using different inverse methods.20,29,30

FIG. 2.

Temperature-density projection of the phase diagram of the (a) diamond forming potential (φdia) and (b) simple cubic forming potential (φsc). Inset (c) zooms in the region around the βSn phase. Symbols and curves are as described in Figure 1.

FIG. 2.

Temperature-density projection of the phase diagram of the (a) diamond forming potential (φdia) and (b) simple cubic forming potential (φsc). Inset (c) zooms in the region around the βSn phase. Symbols and curves are as described in Figure 1.

Close modal

In terms of the qualitative features of the phase diagrams in Figures 1 and 2, we find that there are clear similarities between the two models, which is expected given the corresponding similarities in the pair interactions. Both models feature reentrant melting transitions and many of the same phases including FCC, BCC, low- and high-density hexagonal, and A20 crystals. The fluid-FCC-BCC transitions appearing at low particle densities are similar to those seen in other model soft materials like Hertzian spheres,27 the Gaussian core model,45 and a model for ionic microgels.46 The most striking difference in the phase diagrams of the two models is the presence of the associated target phase—either diamond or simple cubic. Thus, at least for the short-range, soft repulsive potentials considered here, it appears that the potential optimization has a relatively localized effect, having a significant qualitative impact on the phases that appear at intermediate densities.

To summarize, we have used advanced free energy based MC simulations to precisely compute the phase diagrams for two potential models φdia and φsc, which have previously been designed28 to exhibit stable diamond or simple cubic lattice ground states over a wide range of densities. The approach traces out crystal-crystal and fluid-crystal saturation curves within one molecular simulation rather than simulating pairs of phases at several isotherms or isobars to determine the coexistence behavior. We find that the target phases in the models show good thermal stability relative to other potentially competitive crystalline forms, suggesting that our ground-state optimization may also have contributed to the thermal stability of the simple cubic and diamond phases. Despite the simplicity of the model, the overall phase behavior is rich, exhibiting reentrant melting transitions, non-Bravais A20, A7 and βSn crystals, as well as Bravais FCC, BCC, and hexagonal crystalline phases.

Although qualitative features of the optimized models studied here are motivated by the form of effective repulsive forces observed in colloidal and nanoparticle systems, a critical test will be to carry out similar computations on potential models which exhibit a more direct relationship between the model variables and experimentally controllable parameters such as nano-(colloidal) particle type, ligand length, and solvent density. We are currently investigating such a system with experimental collaborators, and we will report on our findings in a future publication.

T.M.T. acknowledges support of the Welch Foundation (F-1696) and the National Science Foundation (CBET-1065357). J.R.E. acknowledges support of the National Science Foundation (CHE-1012356). We also acknowledge the Texas Advanced Computing Center (TACC) at the University of Texas at Austin and the Center for Computational Research at the University at Buffalo for providing HPC resources that have contributed to the research results reported within this paper.

1.
H. S.
Sözüer
and
J. W.
Haus
,
J. Opt. Soc. Am. B
10
,
296
(
1993
).
2.
K. M.
Ho
,
C. T.
Chan
, and
C. M.
Soukoulis
,
Phys. Rev. Lett.
65
,
3152
(
1990
).
3.
C. M.
Soukoulis
and
M.
Wegener
,
Nat. Photonics
5
,
523
(
2011
).
4.
K.
Busch
,
G.
von Freymann
,
S.
Linden
,
S.
Mingaleev
,
L.
Tkeshelashvili
, and
M.
Wegener
,
Phys. Rep.
444
,
101
(
2007
).
5.
C. K.
Ullal
,
M.
Maldovan
,
E. L.
Thomas
,
G.
Chen
,
Y.-J.
Han
, and
S.
Yang
,
Appl. Phys. Lett.
84
,
5434
(
2004
).
6.
S.
Glotzer
and
M.
Solomon
,
Nature Mater.
6
,
557
(
2007
).
7.
O.
Gang
and
Y.
Zhang
,
ACS Nano
5
,
8459
(
2011
).
8.
U.
Agarwal
and
F. A.
Escobedo
,
Nature Mater.
10
,
230
(
2011
).
9.
P. F.
Damasceno
,
M.
Engel
, and
S. C.
Glotzer
,
Science
337
,
453
(
2012
).
10.
C.
Vega
and
P. A.
Monson
,
J. Chem. Phys.
109
,
9938
(
1998
).
11.
H.
Xiong
,
M. Y.
Sfeir
, and
O.
Gang
,
Nano Lett.
10
,
4456
(
2010
).
12.
E. G.
Noya
,
C.
Vega
,
J. P. K.
Doye
, and
A. A.
Louis
,
J. Chem. Phys.
132
,
234511
(
2010
).
13.
F.
Romano
and
F.
Sciortino
,
Nat. Commun.
3
,
975
(
2012
).
14.
G.
Doppelbauer
,
E. G.
Noya
,
E.
Bianchi
, and
G.
Kahl
,
Soft Matter
8
,
7768
(
2012
).
15.
E.
Jankowski
and
S. C.
Glotzer
,
Soft Matter
8
,
2852
(
2012
).
16.
Y.
Wang
,
Y.
Wang
,
D. R.
Breed
,
V. N.
Manoharan
,
L.
Feng
,
A. D.
Hollingsworth
,
M.
Weck
, and
D. J.
Pine
,
Nature (London)
491
,
51
(
2012
).
17.
T.
Vissers
,
Z.
Preisler
,
F.
Smallenburg
,
M.
Dijkstra
, and
F.
Sciortino
,
J. Chem. Phys.
138
,
164505
(
2013
).
18.
S.
Torquato
,
Soft Matter
5
,
1157
(
2009
).
19.
M. C.
Rechtsman
,
F. H.
Stillinger
, and
S.
Torquato
,
Phys. Rev. E
74
,
021404
(
2006
).
20.
M. C.
Rechtsman
,
F. H.
Stillinger
, and
S.
Torquato
,
Phys. Rev. E
75
,
031403
(
2007
).
21.
T.
Yoshida
and
S.
Kamakura
,
Prog. Theor. Phys.
47
,
1801
(
1972
).
22.
F.
Saija
,
S.
Prestipino
, and
G.
Malescio
,
Phys. Rev. E
80
,
031502
(
2009
).
23.
M.
Watzlawek
,
C. N.
Likos
, and
H.
Löwen
,
Phys. Rev. Lett.
82
,
5289
(
1999
).
24.
D.
Gottwald
,
G.
Kahl
, and
C. N.
Likos
,
J. Chem. Phys.
122
,
204503
(
2005
).
25.
T.
Tückmantel
,
F.
Lo Verso
, and
C. N.
Likos
,
Mol. Phys.
107
,
523
(
2009
).
26.
Y. D.
Fomin
,
N. V.
Gribova
,
V. N.
Ryzhov
,
S. M.
Stishov
, and
D.
Frenkel
,
J. Chem. Phys.
129
,
064512
(
2008
).
27.
J. C.
Pàmies
,
A.
Cacciuto
, and
D.
Frenkel
,
J. Chem. Phys.
131
,
044514
(
2009
).
28.
A.
Jain
,
J. R.
Errington
, and
T. M.
Truskett
,
Soft Matter
9
,
3866
(
2013
).
29.
É.
Marcotte
,
F. H.
Stillinger
, and
S.
Torquato
,
J. Chem. Phys.
138
,
061101
(
2013
).
30.
E.
Edlund
,
O.
Lindgren
, and
M. N.
Jacobi
,
J. Chem. Phys.
139
,
024107
(
2013
).
31.
P. S.
Shah
,
J. D.
Holmes
,
K. P.
Johnston
, and
B. A.
Korgel
,
J. Phys. Chem. B
106
,
2545
(
2002
).
32.
P. S.
Shah
,
S.
Husain
,
K. P.
Johnston
, and
B. A.
Korgel
,
J. Phys. Chem. B
106
,
12178
(
2002
).
33.
P.
Schapotschnikow
,
R.
Pool
, and
T. J. H.
Vlugt
,
Nano Lett.
8
,
2930
(
2008
).
34.
L.
Hong
,
A.
Cacciuto
,
E.
Luijten
, and
S.
Granick
,
Nano Lett.
6
,
2510
(
2006
).
35.
É.
Marcotte
,
F. H.
Stillinger
, and
S.
Torquato
,
J. Chem. Phys.
134
,
164105
(
2011
).
36.
The optimized pair potential V(x) as a function of nondimensional interparticle separation x = r/σ has the form
$V(x)= \epsilon \lbrace A x^{-n}\break +\sum _{j=1}^2\lambda _{i}\left(1-\tanh \left[k_{j}\left(x-\delta _j\right)\right]\right) + f_{\text{shift}}(x)\rbrace H[x_{\text{cut}}-x]$
V(x)=ε{Axn+j=12λi1tanhkjxδj+fshift(x)}H[xcutx]
. Here, ε and σ are characteristic energy and length scales chosen so that V(1)/ε = 1,
$x_{\protect \text{cut}}=2.25$
xcut=2.25
is the range, H is the Heaviside step function,
$f_{\protect \text{shift}}(x)\break =Bx^2+Cx+D$
fshift(x)=Bx2+Cx+D
, and constants {B, C, D} are chosen so that
$V(x_{\protect \text{cut}})\break =V^{\prime }(x_{\protect \text{cut}})=V^{\prime \prime }(x_{\protect \text{cut}})=0$
V(xcut)=V(xcut)=V(xcut)=0
.
37.
For the diamond forming potential
$\varphi _{\protect \text{dia}}$
φdia
, the optimized parameters are A = 0.340010893, n = 3.39499,
$\lambda _{\protect \text{1}}=0.732696686$
λ1=0.732696686
,
$k_{\protect \text{1}}=54.1266$
k1=54.1266
,
$\delta _{\protect \text{1}}=2.77156$
δ1=2.77156
,
$\lambda _{\protect \text{2}}=0.605959303$
λ2=0.605959303
,
$k_{\protect \text{2}}=3.71791$
k2=3.71791
,
$\delta _{\protect \text{2}}=1.08107$
δ2=1.08107
, B = −0.0375558, C = 0.20321, and D = 1.17642.
38.
For the simple cubic forming potential
$\varphi _{\protect \text{sc}}$
φsc
, the optimized parameters are A = 0.394620116, n = 5.31842,
$\lambda _{\protect \text{1}}=0.247527755$
λ1=0.247527755
,
$k_{\protect \text{1}} =58.1066$
k1=58.1066
,
$\delta _{\protect \text{1}}=2.63593$
δ1=2.63593
,
$\lambda _{\protect \text{2}}=0.533048319$
λ2=0.533048319
,
$k_{\protect \text{2}}=4.35419$
k2=4.35419
,
$\delta _{\protect \text{2}}=1.05393$
δ2=1.05393
, B = −0.0187544, C = 0.0971676, and D = 0.366054.
39.
D.
Frenkel
and
A. J. C.
Ladd
,
J. Chem. Phys.
81
,
3188
(
1984
).
40.
A. P.
Lyubartsev
,
A. A.
Martsinovski
,
S. V.
Shevkunov
, and
P. N.
Vorontsov-Velyaminov
,
J. Chem. Phys.
96
,
1776
(
1992
).
41.
K. S.
Rane
,
S.
Murali
, and
J. R.
Errington
,
J. Chem. Theory Comput.
9
,
2552
(
2013
).
42.
J. R.
Errington
,
J. Chem. Phys.
118
,
9915
(
2003
).
43.
See supplementary material at http://dx.doi.org/10.1063/1.4825173 for simulation results of assembly to the corresponding target structures upon isochoric quenching.
44.
S.
Prestipino
,
F.
Saija
, and
G.
Malescio
,
Soft Matter
5
,
2795
(
2009
).
45.
S.
Prestipino
,
F.
Saija
, and
P. V.
Giaquinta
,
Phys. Rev. E
71
,
050102
(R) (
2005
).
46.
D.
Gottwald
,
C. N.
Likos
,
G.
Kahl
, and
H.
Löwen
,
Phys. Rev. Lett.
92
,
068301
(
2004
).

Supplementary Material