Cluster Crystals Stabilized by Hydrophobic and Electrostatic

Article Cite This: J. Phys. Chem. B 2018, 122, 2669−2682

pubs.acs.org/JPCB

Cluster Crystals Stabilized by Hydrophobic and Electrostatic Interactions A. Baumketner,*,† A. Stelmakh,‡ and W. Cai§ †

Institute for Condensed Matter Physics, NAS of Ukraine, 1 Svientsistsky Street, Lviv 79011, Ukraine Department of Chemistry, Ivan Franko Lviv National University, 6 Kyrylo and Mefodii Street, Lviv 79005, Ukraine § Department of Mathematics, Southern Methodist University, Dallas, Texas 75252, United States

Downloaded via KAOHSIUNG MEDICAL UNIV on June 30, 2018 at 17:54:56 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

‡

ABSTRACT: Cluster crystals are crystalline materials in which each site is occupied by multiple identical particles, atoms, colloids, or polymers. There are two classes of systems that make cluster crystals. One is composed of particles that interact via potentials that are bound at the origin and thus are able to penetrate each other. The other consists of noninterpenetrating particles whose interaction potential diverges at the origin. The goal of this work is to find which systems of the second class can make cluster crystals that are stable at room temperature. First, the general properties of the required potentials are established using an analytical model and Monte Carlo simulations. Next, we ask how such potentials can be constructed by combining hydrophobic attraction and electrostatic repulsion. A colloid model with a hard-sphere core and a repulsive wall is introduced to mimic the hydrophobic interaction. Charge is added to create long-range repulsion. A search in the parameter space of the colloid size, counterion type, and charge configuration uncovers several models for which effective colloid−colloid interaction, determined in explicit solvent as a potential of mean force, has the necessary shape. For the effective potential, cluster crystals are confirmed as low free-energy configurations in replica-exchange molecular dynamics simulations, which also generate the respective transition temperature. The model that exhibits a transition above room temperature is further studied in explicit solvent. Simulations on a 10 ns time scale show that crystalline conformations are stable below the target temperature but disintegrate rapidly above it, supporting the idea that hydrophobic and electrostatic interactions are sufficient to induce an assembly of cluster crystals. Finally, we discuss which physical systems are good candidates for experimental observations of cluster crystals.

1. INTRODUCTION The term “cluster crystal” (CC) refers to crystalline materials in which each site is occupied by multiple entities of the same kind, atoms, colloids, or polymers that, collectively, are known as clusters. Because of their exciting novel properties, such systems are of great interest to fundamental research. At the same time, strong periodic order makes them a suitable candidate for use in various technological applications, for instance, nanopatterning.1 Much of what we know today about cluster crystals comes from theory and computer simulations. Historically, the first observation of a CC was made in a Monte Carlo (MC) simulation2 of the model system interacting via a step potential with a finite value at the origin. The studied potential is repulsive, but allows for a full overlap between particles. In a subsequent work (see ref 3 and references therein), these two features were seen as central to the ability of the system to assemble into a large array of various mesoscopic structures, including cluster crystals, columnar crystals, or lamellae. The assembly process is driven by the need of each particle to minimize the number of its immediate neighbors. It turns out that this can be accomplished more easily in crystalline ordered structures than in the disordered glasslike ones. A quantitative criterion2,4 based on the Fourier transform of the given © 2018 American Chemical Society

repulsive potential can be used to judge whether or not it favors ordered structures. Over the last 10 years, a large number of studies have focused on cluster crystals in repulsive systems.5−7 Some of the aspects addressed in these studies include: (a) development of new free-energy methods,8 (b) collective dynamics,9 (c) diffusion processes,10 and (d) behavior under shear.11 A major effort was also made to test the predictions made with the help of the model potentials in real-life systems. Great strides in this regard were made for polymer systems,12,13 which are known to allow overlap of constituent monomeric chains. Dendrimers have been of particular interest due to their accessibility in both experimental and simulation settings. A recent tour de force simulation14 reported a dendrimer model with monomer resolution spontaneously assembling into a cluster crystal. This is a major step forward, potentially paving the way for experimental observation of these systems, which so far has remained elusive.13 On the other hand, repulsive overlapping particles are not the only ones capable of making cluster crystals. It was briefly noted in previous studies 2 that step potentials with Received: November 27, 2017 Revised: February 9, 2018 Published: February 12, 2018 2669

DOI: 10.1021/acs.jpcb.7b11662 J. Phys. Chem. B 2018, 122, 2669−2682

Article

The Journal of Physical Chemistry B

distances and a repulsive shoulder at long distances. We next ask whether it is possible to design colloidal particles so that they interact via effective potentials closely matching those that have been identified in the model systems. For this purpose, we employ a scheme in which hydrophobic interactions are chosen to represent the short-range attraction, whereas electrostatic interactions are used to create the long-range repulsion. The colloids are modeled as hard particles of certain radius with repulsive walls. Effective potentials are extracted as potentials of mean force (PMF) in explicit solvent simulations. With no charge, the colloids exhibit an attractive well at short distances, demonstrating their hydrophobic character. The addition of charge has two effects: (a) it creates a repulsive tail and (b) it shifts the position and the value of the minimum. Changing the magnitude and location of the charge within the colloid allows us to manipulate the shape of the resulting effective potential. An exhaustive search yields several models, in which the charge is either concentrated in one point inside the colloid or distributed on its surface, that generate PMFs similar to those identified for the generic model. These potentials are then tested in replica-exchange simulations to confirm that they lead to cluster crystals as the most stable phase at room temperature. Further analysis is conducted only for those potentials that generate small clusters, 13 colloids or less. The salient features of such systems is a minimum of ∼10 kJ/mol and a maximum of ∼5 kJ/mol, whereas the radius of the colloid varies between 3 and 7 Å. As the description of the system in terms of effective potentials, or implicit solvent, is approximate, further tests are conducted in more accurate explicit solvent simulations. For the selected few models, cluster crystals are shown to be stable at room temperature on the time scale of 10 ns. An increase in temperature leads to the crystal meltdown with the transition temperature estimated from explicit solvent model to be slightly above that predicted in the effective-potential approach. The error is remarkably small, within 10%, showing that the nonadditivity effects are not strong for the studied systems. Our study is the first to demonstrate that a combination of hydrophobic and electrostatic interactions is sufficient to support cluster crystal in appropriately chosen systems. We provide specific recipes for how to place charges on a colloidal particle for it to be able to assemble into cluster crystals, which should make the search for experimentally observable systems much easier.

impenetrable hard cores could also produce these structures. More recent studies of the same15,16 or similar finite-size repulsive system17 confirm these conclusions. Nonoverlapping particles interacting via short-range attractive and long-range repulsive (SALR) potentials18 are also able to make cluster crystals. The attractive part in the potential induces a phase transition into the liquid state.19 The repulsive part frustrates that transition by limiting the size to which a nascent liquid nucleus can grow. For sufficiently strong repulsions, the liquid phase is preempted by a number of finite-size structures known as microphases. A large number of microphases have been identified, which can be both ordered and disordered, with cluster crystals being one of them. The first details of how microphase separation takes place were gleaned in early theoretical studies.20−22 The main difference with the purely repulsive potentials is that SALR clusters are energetically favorable. This has consequences for the cluster size distribution as well as for the scaling of the cluster size with the number of particles: SALR clusters are compact and sharply peaked around certain size (which is density-dependent) in contrast to the clusters in repulsive systems, which are extended and come in a broad variety of sizes.23 At sufficiently high temperatures, preassembled clusters in SALR systems remain in equilibrium with monomers, forming a cluster fluid. The fluid may experience freezing into a number of microphases depending on the system density as the temperature is lowered. With increasing density, the sequence of available phases appears to be universal24 and begins with crystals before continuing on to columns and to lamellae. If cooling is conducted too fast, kinetic arrest may prevent the appearance of an ordered state favoring instead glassy cluster states. Theoretical predictions are available for the transition temperature into both a repulsive glass and a cluster crystal state.25 The boundaries of the ordered microphases can also be delineated by the so-called λ-line, determined within the meanfield approximation.26 In-depth characterization of the microphases populated by a particular potential remains a huge challenge. Theoretical treatments based on continuum approximation or integral equations are efficient, but may lack the necessary level of accuracy. Computer simulations, on the other hand, are accurate but also computationally expensive. Adding to the problem is the slow dynamics that accompanies phase transitions. A combination of these factors has contributed to the fact that out of a limited simulation work on SALR systems27−39 only a few papers32,35−39 were actually able to observe ordered microphases directly. Most of the studies focused on columns and lamellae,32,35−39 whereas cluster crystals were seen only in two recent reports.37,38 To the best of our knowledge, no experimental observations of cluster crystals have been made for such systems to date. At present, cluster crystals in systems represented by potentials of overlapping particles are much better studied than their nonoverlapping counterparts, a gap that this paper aims to fill. We consider particles that are not allowed to interpenetrate and ask what it would take for them to make a stable cluster crystal at room temperature. No preassumptions are made regarding the shape of the interaction potential, SALR or otherwise. Instead, a thorough search by computer simulations is conducted in the parameter space of a generic interparticle potential introduced by us earlier.23 A number of potentials with the desired property are uncovered, all of which display either a local or a global attractive minimum at short

2. METHODS 2.1. Analytical Potential. We used analytical potential introduced in our prior work23 to investigate conditions under which cluster crystals can be observed. The potential has a hard wall, a triangular minimum describing attraction, and a Yukawa tail describing repulsion between particles. The exact formula u(r ) = ⎧+∞ , r < σH ⎪ ⎪−ar + b + ε , σ < r < σH + σA , a = 2δε w H ⎪ 2 σA − σH ⎪ σH + σA σ + σH ⎨ < r < σA , b = δε A ⎪ ar − b + εw , 2 σA − σH ⎪ ⎪ −κ(r − σA ) e ⎪ , σA < r ⎪ εbσA ⎩ r

(1)

depends on six parameters, σH, σA, εw, εb, δε, and κ. Specific values are discussed in the main text for each particular system. 2670

DOI: 10.1021/acs.jpcb.7b11662 J. Phys. Chem. B 2018, 122, 2669−2682

Article

The Journal of Physical Chemistry B

Figure 1. Details of hydrophobic interactions between hard-sphere colloids with repulsive walls. (a) PMF computed for varying radius of the colloid R. For all sizes, the potential turns zero at approximately r = 2R + 8 Å. The inset shows the hydrophobic solvation cost, solid symbols, as a function of R. The curves denote different approximations discussed in the text. (b) Different interpretations of the solvent-accessible surface area (SASA). The blue curve denotes area measured on the surface of the colloid, and the red curve shows the surface swept by the center of a probe molecule.

In all systems, δε = 1.9 kJ/mol. To perform a search in the parameter space, we used Monte Carlo (MC)-simulated annealing runs. A system with N = 216 particles was considered. We tested larger systems, but the results differed little other than in the fact that obtaining a crystalline structure becomes very difficult for N > 512. Temperature was lowered from 700 to 50 K over 10 different values constructed using a constant increment in the inverse temperature. At each temperature, the systems were simulated for 2 × 106 MC sweeps. For each system, at least 10 repeat runs were performed. For systems that give rise to cluster crystals, at least 100 repeats were done. The density varied for each system so as to cover the area where clusters of size 3 < Nc < 15 are observed. Parameters εw and εb were changed in 1 kJ/mol increments. The inverse screening length κ changed in 1 nm−1 increments. Two main values, 2σH and 5σH, were used for σA. Multiple other values were investigated for specific sets of parameters, in particular, σA = 1.8σH and σA = 2.2σH. For those parameter sets that lead to cluster crystals, additional simulations by the replica-exchange40−43 algorithm were performed. A total of 28 replicas were considered, spaced exponentially from 220 to 500 K. Simulations were run for 30 × 106 MC steps. The acceptance probabilities between neighboring replicas were 10% or more. The generated trajectories were then subjected to multiple-histogram analysis44,45 to extract relevant thermodynamic properties as a function of temperature. Simulations were started from crystalline initial structures. Tests performed on different initial structures showed that structural ensembles populated at temperatures below Tc are not converged. For instance, we were unable to determine the population of different crystalline structures in this way. It seems that this method converges too slowly to be suitable for this task. 2.2. Hard-Sphere Colloid with Repulsive Walls. Hydrophobic particles in water have the tendency to stick together. We derived a model that properly captures this effect, while remaining computationally affordable and not specific to any particular system. The starting point of the derivation is a hardsphere particle of radius R that interacts with any other point particle in the system through a repulsive potential. Assuming

12

( σr )

that the Lennard-Jones repulsion, 4ϵ

, is used and the

surface of the colloid is uniformly covered in interaction points, one can arrive, after integration, at the following formula for the solute−point (SP) interaction uSP(r ) =

12 ⎞ 4ϵSPσSP 1⎛ 1 1 , − ⎜ 10 10 ⎟ 20 Rr ⎝ (r − R ) (r + R ) ⎠

⎛ σ ⎞12 lim uSP(r ) = 4ϵSP⎜ SP ⎟ ⎝ r ⎠ R→0

(2)

where r is the distance from the center of the colloid to the point. By design, the colloid shrinks to an interaction point when parameter R tends to zero, so the potential recovers the proper form in that limit. For finite R, it describes interactions of colloids with water, uSW(r), or counterions, uSI(r). The same principle of spreading repulsive points on the surface is used to obtain the following expression for the solute−solute (SS) interaction uSS(r ) =

12 4ϵSSσSS 1 (65536R16 − 147456R14r 2 5 9r10(r 2 − 4R2)9

+ 147456R12r 4 − 86016R10r 6 + 32256R8r 8 − 7488R6r10 ⎛ σ ⎞12 + 2688R4r12 + 360R2r14 + 45r16), lim uSS(r ) = 4ϵSS⎜ SS ⎟ ⎝ r ⎠ R→0

(3)

where r is the distance between the centers of the colloids. Here, again the interaction between two point particles is recovered in the limit of R → 0. The relevant parameters, σSS = 0.35 nm and ϵSS = 0.276144 kJ/mol, were adapted from the allatom optimized potentials for liquid simulations (OPLS-AA) force field,46 alkane carbon group, so in the limit of R → 0, formula 3 describes repulsion between two methane molecules. Geometric combination rule47 is applied to generate van der Waals cross terms. For instance, one obtains σSW = 0.33207 nm and ϵSW = 0.4192 kJ/mol for colloid−water interactions, where TIP3P48 parameters were used for the water. A similar calculation was performed for the colloid−counterions potentials. 2671

DOI: 10.1021/acs.jpcb.7b11662 J. Phys. Chem. B 2018, 122, 2669−2682

Article

The Journal of Physical Chemistry B 2.3. Potential of Mean Force Computation. The effective interactions between colloids u(r), which model the solvent implicitly, were computed by the umbrella sampling method.49 The colloids of radius R were solvated in a box of water with dimensions 2R + 20 Å in y and z directions and 2R + 40 Å in x direction. Umbrella potential was applied along x axis with the force constant k = 2000 kJ/mol nm2. A total of 15 windows were considered, spaced at locations between 2R + 1 Å and 2R + 16 Å with the increment of 1 Å. Simulations were run for 3 ns for each window. Weighted histogram analysis method was used to obtain unbiased distribution functions.44,50 The effective interactions are shown in Figure 1a for varying R. To make comparison easier for different systems, the data are plotted as a function of the distance between colloids d = r − 2R. All curves have the same shape with a single minimum and no other extrema. The position of the minimum varies slowly with R, decreasing from d = 2.3 Å for R = 1 Å to d = 1.0 Å for R = 7 Å. Interestingly, the point at which potentials turn zero seems to be the same for all systems, 2R + 8 Å, which corresponds to the 8 Å solvation layer between the colloids. The cost of the hydrophobic attraction between two colloids, UHB(R), can be determined from the minimum value of the potential. It is seen from the insets of Figure 1a that UHB(R) grows rapidly with R. A power law dependence 5.9 + 2.8x1.55, shown by the red curve in the figure, provides a good fit to the data. A plot on logarithmic scale in the inset provides additional details to help judge the quality of this approximation. For the sake of comparison, we also analyzed our data in reference to a popular model of hydrophobic solvation based on solvent-accessible surface area (SASA).51 We note that there is some ambiguity in the definition of SASA, illustrated in Figure 1b for the case of two colloids with radius R and separated one from another by distance 2R + d. If a solvent molecule has a finite radius w, geometric constraints will prevent it from entering the interstitial space between the colloids if d < 2w. As a result, certain excluded, inaccessible to solvent, region on the colloid will be formed, defined by a corresponding solid angle. As illustrated in Figure 1b, the surface area of the excluded zone, shown by the blue curve, is

(

ΔS1 = 2πR2 1 −

R + 0.5d R+w

the more realistic soft-particle model is adopted, w includes both the radius of the solvent and the radius of the interaction points located on the colloid surface. The latter can be taken as 0.5σSW = 1.7 Å, which leads to a more reasonable estimate of the water radius Rw = 1.9 Å. We note that the maximum in the colloid−solvent pair distribution function is at a distance of 3.1 Å away from the colloid surface. So, the fitted radius seems to point to the first solvation shell of the colloid, supporting the loss of solvent in that shell as the main mechanism underlying hydrophobic attraction. An alternative definition of the inaccessible surface is to take the surface swept by the center of a probe molecule, as illustrated in Figure 1b by ΔS2 and the red curve. The hydrophobic solvation then reads ⎛ R + 0.5d(R ) ⎞ ⎟ ΔUHB(R ) = γ 2ΔS2 = γ 4π(R + w)2 ⎜1 − ⎝ R+w ⎠

where γ has the same meaning as before. Fitting produces γ = 666 J/mol/Å2 and w = 1.3 Å. The corresponding UHB(R) is shown in Figure 1a by the blue curve. The fit is notably worse than for the other two models, but overall not bad. The surface tension is too high compared to the literature, whereas the probe radius is too small. For the common choice of Rw = 1.4 Å,53 for instance, the expected w = 3.1 Å. Furthermore, it is not clear how to interpret the surface delineated by w = 1.3 Å. As no solvent molecules are allowed to approach the colloid at such short distances, the surface traced out by the center of the solvent molecules is clearly not suitable. It appears, therefore, that although both SASA definitions can be parameterized to reproduce the simulation data well, the model relying on the surface area of the actual colloid has fewer inconsistencies. To compute PMFs for colloids carrying charge, we used the same setup as described above. The van der Waals parameters σ = 0.325 nm and ϵ = 0.71128 kJ/mol for the charges on the colloid surface were borrowed from the peptide nitrogen atom in OPLS-AA force field.46 Counterions were added to neutralize the total charge of the simulation box. Because counterion concentration c has a strong effect on the resulting effective potential, care must be taken to ensure that the same c is used in both PMF simulations in which u(r) is computed, and implicit solvent simulations using u(r), in which cluster crystals are studied. Obviously, this cannot be accomplished exactly, as the concentration in the cluster crystal is not known ahead of time. Fortunately, the concentration dependence of the effective potential can be inferred from an analytical model for long distances between colloids and low salt concentration72

). The loss of solvent incurs a loss of

internal energy due to the solvent−colloid interactions, which, in the present model, is proportional to ΔS1. Assuming that the free energy can be approximated by the internal energy (entropic effects are neglected), hydrophobic solvation can be written as ΔUHB(R ) = UHB(R ) − UHB(0) = γ 2ΔS1 ⎛ R + 0.5d(R ) ⎞ ⎟ = γ 4πR2⎜1 − ⎝ R+w ⎠

u(r ) = f

2 qeff

r

e−κ(r − rD)

(4)

⎡ kJ Å ⎤ where f = 138.935/8⎣ mol e2 ⎦, r is measured in Angstroms, qeff is the effective charge of the colloid, which includes the real and the screening charge, κ is the inverse screening length, and rD is a certain cutoff distance from which a screened exponential decay in the PMF begins. According to the Debye formula,72 the screening constant κ is related to ni, the counterion

where γ is a proportionality constant, or surface tension, and d is assumed to be explicitly dependent on R. Of course, this approximation can be accurate only in the limit of large R. It contains two undefined parameters γ and w, which can be determined by fitting. If we take UHB(0) = 5.9 kJ/mol, then γ = 318 J/mol Å2 and w = 3.6 Å for our data and the corresponding UHB(R) is shown in the inset of Figure 1a by the green curve. The surface tension agrees well with the value of 301 J/mol/Å2 reported in the literature,52 whereas UHB(R) is accurate for large R > 3 Å, as intended. The fitted radius 3.6 Å seems unreasonably large, however. This is the consequence of representing solvent molecules as hard-sphere particles. When

concentration, as κ =

2 4πf n iqi 10 kT

, where k is the Boltzmann

constant, T is the temperature, and qi is the charge of counterions. Therefore, if u(r) is available for one concentration, it can be approximated for another concentration by 2672

DOI: 10.1021/acs.jpcb.7b11662 J. Phys. Chem. B 2018, 122, 2669−2682

Article

The Journal of Physical Chemistry B

Figure 2. Results obtained with the help of a model potential. (a) Schematics of potentials leading to cluster crystals. (b−f) Data for the model termed type (c), shown by the red line in (a). (b) Average cluster size Nc as a function of particle packing fraction η. (c) Plot of distribution over cluster size s for varying temperature and packing fraction η = 0.034. (d, e) Data for this η. (d) The average value of bond-orientational order (BOO) parameter Q6, and its dispersion ⟨ΔQ26⟩, as a function of temperature. Temperature of the transition into the cluster crystal state is indicated by the broken line. (e) A snapshot from the trajectory at T = 272 K, demonstrating the degree of disorder in the structure. (f) Transition temperature as a function of cluster size Nc. Clusters of size 8 and larger lead to crystals that are stable at T = 300 K.

molecule to −1. Parameters for propionate ions were constructed in a similar manner. 2.6. Details of All-Atom Crystal Simulations. Initial crystal configurations were constructed manually using the colloid density and cluster size suggested by the implicit solvent simulations. Electrostatic interactions were computed by the particle-mesh Ewald method.58 TIP3P48 model of water was used for solvent. Water geometry was kept fixed by the SETTLE59 algorithm. Bond lengths among colloid charges were kept fixed by the LINCS60 algorithm. The time step in all simulations was 2 fs. Constant temperature was maintained by the Nosé−Hoover thermostat61 with a 0.5 ps time constant. The duration of the simulations varied with temperature but always exceeded 10 ns.

adjusting the tail of the potential. The other parameters needed in this procedure, qeff and rD, can be obtained by fitting and are assumed to be little affected by small concentration changes. 2.4. Simulations Using Effective Potentials. Systems interacting via PMF potentials were studied using the same protocol as those of the analytical potential. First, simulated annealing runs by the molecular dynamics (MD) method were conducted to find the low free-energy conformations and to determine the dependence of the mean cluster size Nc on density. Later on, for those systems that make cluster crystals, additional replica-exchange simulations were conducted. A total of 16 replicas were considered at temperatures spaced between 200 and 450 K. The simulations were run for 10 ns and yielded acceptance probabilities between replicas not less than 10%. Conformations for all replicas were clustered for the computation of Q6. Multiple-histogram analysis44,45 was applied to compute Q6 and its dispersion as a function of temperature. The maximum in the latter quantity was used to identify transition temperature into the crystalline state. The GROMACS software package54 was used to conduct all MD simulations. 2.5. Counterions. We tested a number of counterions in the present work. Both positive and negative charges were considered. Parameters from Jorgensen’s group were used for halide F− and Cl− ions55 and for I−.5656 Parameters from Ȧ qvist’s work57 were adopted for alkali Na+, Li+, and K+ ions. The OPLS-AA force field46 was used to model acetate ions. Specifically, geometry and parameters of the bond-angle, dihedral-angle, and improper-angle potentials of glutamic acid were employed. Also all charges were copied from this amino acid residue except for the alkane carbon whose charge was changed from −0.22 to −0.28 to bring the total charge of the

3. RESULTS 3.1. Potentials Leading to Cluster Crystals. We used a generic-type potential introduced in our prior work23 to characterize conditions under which cluster crystals can be observed. The potential represents a hard-sphere particle of diameter σH that has a triangular minimum with the lowest point of εw and extending from σH to σA. The minimum is followed by an energy maximum of εb, after which there is a decline to zero according to the Yukawa, or the screened 1 electrostatics law, ∼ r e−κ(r − σA), where κ is the inverse screening length. The five parameters, σH, σA, εw, εb, and κ encode a whole plethora of potential shapes, including pure repulsive potentials and SALR potentials.23 These parameters were varied systematically in Monte Carlo simulations to determine which of them favor cluster crystals. Because of the large scope of the required search, a two-stage approach was adopted. First, 2673

DOI: 10.1021/acs.jpcb.7b11662 J. Phys. Chem. B 2018, 122, 2669−2682

Article

The Journal of Physical Chemistry B simulated annealing runs were conducted to find possible candidates for stable low-temperature structures. Structures with crystalline order were selected for further analysis. In the second stage, the structures that displayed an apparent transition temperature above 300 K were studied by more accurate replica-exchange simulations. More details are reported in the Methods section. The selected potentials were screened against two criteria. First, they needed to display a transition temperature above 290 K so that crystals are stable at room temperature. Second, the clusters needed to be small. This second requirement is motivated by the fact that only in small clusters all particles are located on the surface, whereas in large clusters both bulk and surface particles are present. Although interactions among surface particles can be reasonably well approximated within the effective pair potential approach, the method that is used in this work to compute intercolloid potentials, particles buried in the interior of a large object are subject to a substantial error due to nonadditivity effects.62 To avoid such effects, we limited the size of the considered clusters to 13 particles, the largest size before bulk particles begin to emerge according to our visual inspections. This estimate is well within the range reported in the literature for atomic clusters.63,64 All potentials that meet the above two criteria can be divided into three groups, as explained in Figure 2a. The first group, labeled type (a), consists of potentials with a shallow local minimum and a high energy barrier. The width of the minimum has to be more than twice the size of the hard core diameter, σA = >2.2σH, whereas the height of the energy barrier has to be at least 14 kJ/mol. A local minimum is also present in the second group, type (b). It is now allowed to be much narrower, σA = >1.8σH, but must exhibit a more significant depth Δε = εb − εw of at least 10 kJ/mol. The third group, type (c), is similar to group (b), but its minimum is global instead of local. The screening length has little effect on the transition temperature in all groups as long as it is in the range defined by the size of the particle, 1/σA < κ < 1/σH. For illustration purposes, we present in Figure 2b−f the results of a complete thermodynamic analysis for the potential belonging to type (c) and defined by the following parameters: σH = 0.3 nm, σA = 0.54 nm, εw = −5 kJ/mol, εb = 8 kJ/mol, and κ = 4 nm−1 (see Figure 2a). Simulated annealing runs uncovered cluster formation for this potential at low temperature and varying density. A particle was considered part of an existing cluster if its distance to any other member of the cluster is less than or equal σA. Figure 2b shows average cluster size Nc obtained in simulated annealing runs at the lowest temperature T = 80 K. Here, the particle density ρ is expressed in terms of the volume fraction η = ρσ3Hπ/6. The curve has a concave shape with a clear linear dependence for η > 0.03, in agreement with earlier reports for SALR potentials.25,65 The minimal-energy conformation identified in our simulations is face-centered cubic (FCC) crystal. Figure 2c shows fraction of particles P(s) involved in clusters of size s obtained in a replica-exchange simulation for η = 0.034 that corresponds to Nc = 8. The initial configuration for this simulation was assembled manually from 256 particles that were put into 32 identical clusters and placed in an FCC arrangement inside a cubic simulation box of length L = 4.76 nm. The temperatures were distributed across 28 replicas according to the geometric progression in the range of 220−500 K. It is seen that at T = 220 K there are three nonzero values in P(s) that correspond to s = 7, 8, and 9. The particles, therefore, begin to migrate across

available crystal sites at a very low temperature. As a consequence, the low-temperature phase is not an ideal crystal with equal number of particles at each site but a crystal with a degree of polydispersity, in which a certain distribution is established in the size of each cluster. As the temperature is increased, the distribution over the cluster size becomes broader. The broadening is most dramatic around T ∼ 300 K, indicated in Figure 2c by red arrows, where a jump in the cluster size from 10 to 29 is observed. Above this temperature, P(s) has a multimodal appearance with maxima located at s = 8, 16, and 24, which is consistent with the system being in a cluster fluid state composed of octamer as the main, or preferred, cluster and its multimers created by the cluster− cluster associations. Cluster multimers survive until T = 400 K, at which point, P(s) acquires a shape with only one broad maximum. To learn more about the structural transformations taking place in the system with temperature, we analyze the behavior of bond-orientational order (BOO) parameters.66 Reporting on the local order around each cluster, BOOs can reliably measure the difference among different types of lattices as well as disordered glassy states. We specifically focus on average Q6 as this parameter seems to work well at distinguishing between crystalline and glassy structures.67 All conformations in our simulations were clustered. Geometrical centers of the generated clusters, or centroids, were then used to compute Q6. Centroids separated by less than 2 nm (a distance slightly longer than 1.7 nm, the position of the maximum in the centroids pair distribution function) were considered nearest neighbors in this computation. Figure 2d shows average values of Q6, as well as its dispersion ⟨ΔQ26⟩, as a function of temperature, computed using the multiple-histogram analysis method.44,45 Starting at a value of 0.47 at T = 200 K, Q6 declines gradually with temperature to a value of 0.32 at T = 450 K. As an ideal FCC crystal has Q6 = 0.57 and values near 0.3 are typically seen for disordered states,67 it is clear that the initial crystal melts with temperature. The melting point appears to be around T = 300 K, where Q6 experiences the fastest decline. More accurately, the transition temperature Tc can be estimated from the position of the maximum in the dispersion curve, plotted in the inset of Figure 2d and showing that Tc = 310 K. The order parameter at this temperature is 0.38, suggesting that this value can be used as the cutoff parameter with which to sort out crystalline and disordered states. Figure 2e shows a conformation observed in the trajectory at T = 272 K, which is slightly below the transition temperature. The crystalline structure is clearly visible although the system is far from an ideal crystal. Multiple deformations and dislocations lead to the order parameter of the reported structure, 0.46, dropping below 0.57 of the ideal crystal, yet still remaining above the crystalline cutoff threshold. The structure contains 22 clusters of size s = 8 and 5 clusters of size 7 and 9 each, proving that the low-temperature phase is polydisperse. The analysis presented in Figure 2c−e was repeated for a number of densities. The resulting transition temperature is shown in Figure 2f as a function of the average cluster size Nc. It can be well approximated by a line for Nc ≥ 8, whereas for smaller clusters, it declines rapidly. It is clear from the figure that for a select target temperature Tg there is a minimal cluster size N*c , corresponding to a minimum particle density η*, above which cluster crystal is stable. For example Nc* = 7.8 for Tg = 300 K, which corresponds to the particle density η* = 0.032. 2674

DOI: 10.1021/acs.jpcb.7b11662 J. Phys. Chem. B 2018, 122, 2669−2682

Article

The Journal of Physical Chemistry B 3.2. Potentials for Hydrophobic Particles Carrying Charge. 3.2.1. General Considerations. 3.2.1.1. Interaction Schemes. In regard to Figure 2a, the next pertinent question to ask is how to realize potentials shown there in actual physical systems. All identified potentials have an attractive and a repulsive part. There is more than one way to achieve this shape. The simplest systems, perhaps, are mixtures of hardsphere particles of differing size. The larger particles experience interaction due to the presence of the smaller particles, the socalled depletion force, which under appropriate conditions exhibits one minimum and one maximum,68 similarly to some of the potentials shown in Figure 2a. As demonstrated both experimentally69 and theoretically,70 the desired shape may also arise as a consequence of the low dimensionality of the physical space. In this work, we pursue a third possibility and use the potential derived previously by our group for protein lysozyme23 as a cue to introduce a scheme in which hydrophobic interactions are employed to model the attractive force, whereas electrostatic interactions are used to represent the repulsive force. In support of this combination, we note that the chosen two interactions are central to key processes in chemistry and biology.71 3.2.1.2. Bounds on the Parameters of the Colloids. Initial estimates for the parameters of colloids making cluster crystals can be obtained from analytical theory. As described in the Methods section, a hard-sphere particle of radius R is introduced as a model of hydrophobic interactions. Charge q is placed at the center of the particle. Let us now assume that we aim to construct with this model an SALR potential that has a global minimum, as shown, for example, in Figure 3. Let us

electron charge, and the dielectric constant is assumed to be 80 here and elsewhere in the paper. There are two conditions for the total energy to meet: (a) the energy at the contact distance has to be negative to make sure that the minimum is global, and (b) the height of the barrier needs to be ∼10 kJ/mol, according to our estimates from Figure 2a. Let us assume that the maximum occurs where the hydrophobic potential turns zero, 2R + 8 Å, whereas the minimum appears at 2R + 2 Å, which is the diameter plus a certain distance allowed between the colloids because of the soft wall. In this case, the two conditions translate into two inequalities,

fq2 2R + 8

> 10 a n d

2

fq 2R + 2

< 6 + 2.8R3/2 , which can be satisfied simultaneously only for R > Rmin = 3 Å. Thus, there exists a lower bound on the size of the colloids. The charge for the allowed radii is found in

the range

10(2R + 8) f

< q2
rD = 1.8 nm can be well fitted by the model (eq 4), in which qeff is treated as an adjustable parameter, whereas a value of 0.82 nm−1, obtained from the Debye formula, is used for κ. The analytical extension is shown in the figure by solid line. Similar extensions were obtained for other potentials discussed in this work. Further increase in q leads to further shifts of the PMF upward, whereas the minimum gradually diminishes. For q = 3e, a local minimum is present, but its depth of 1 kJ/mol is so negligibly small that the entire potential

Figure 3. Hypothetical potential that can be obtained by combining short-range hydrophobic attraction with long-range electrostatic repulsion. The distance at which the colloids are in contact corresponds to the minimum. The maximum is found where the forces due to the two interactions are same in magnitude but opposite in sign, which can be approximated by the distance where the hydrophobic attraction vanishes.

further assume that the full interparticle potential is the algebraic sum of its two components, hydrophobic attraction and electrostatic repulsion (little influence of one on another). The hydrophobic energy released upon making a contact can be well described, according to our model, by the expression uHB(R) = 6 + 2.8R3/2, where R is measured in angstroms and kJ/mol is used for the unit of energy. The electrostatic repulsion at long distances can be approximated by the q2

Coulomb law in a dielectric medium uel(r ) = f r , where ⎡ kJ Å ⎤ f = 138.935/8⎣ , r is measured in angstroms, e is the mol e 2 ⎦ 2675

DOI: 10.1021/acs.jpcb.7b11662 J. Phys. Chem. B 2018, 122, 2669−2682

Article

The Journal of Physical Chemistry B

Figure 4. Potentials of mean force obtained in this work for various systems. (a) Results for the colloid of size R = 5 Å and central placement of the charge, simulated in the presence of F− ions. The maximum charge for which SALR shape is seen is q = 2e. (b) Effect of shifting the charge off center by 2 Å for q = 3e, the same size of the colloid and Cl− counterions. (c) Similar data for a larger colloid, R = 6 Å, and two charges, q = 3e and 4e, placed 2 Å off center.

kJ/mol, which is higher than the value of 2.5 kJ/mol attainable for R = 5 Å. However, together with the barrier, the minimum also goes up, as is clear from the potential for R = 6 Å and q = 3e. The dependence of the potential on R and q is very complex and requires a substantial effort by simulations to identify an appropriate set of parameters. 3.2.2.3. Cluster Crystals in the One-Charge Model. We examined a total of 109 models, of which 60 contained a central charge and the remaining 49 had the charge displaced. Our search successfully identified systems capable of making cluster crystals. Figure 5 shows the PMF of one such system in which

can be considered repulsive. A repulsive and very similar potential is seen for q = 4e, making it evident that qmax = 2e is the maximum charge for which an SALR potential is seen. The maximum charge is found to correlate with the colloid size: in smaller colloids, qmax is smaller and vice versa in larger colloids. For R = 3 Å, qmax equals 1e, whereas for R = 7 Å, it is 3e. For the R = 5 Å case, a comparison for q = 3e and 4e reveals that the two curves begin to converge. The reason behind the trend toward convergence is the excessive binding of counterions to the colloid. When q is large, charging the colloid by one extra unit results in one additional counterion bound to it. As a consequence, the total screened charge remains almost unchanged, causing the PMF to vary very little. The barrier in the SALR potential for q = 2e is too low for cluster crystals, ruling out potentials of type (c) for the R = 5 Å model. The barrier in the potentials with a local minimum, see q = 3e and 4e, is also too low, implying that the type (b) and (a) shapes cannot be attained in this model either. In addition, the minimum becomes notably narrower with the increase of q, suggesting that type (a) shape cannot be reached even for larger R. 3.2.2.2. Off-Center Charge. As an extension of the model, we studied the effect of shifting the charge away from the colloid’s center. In two colloids in a close-contact configuration, displacement of the charge sites increases the distance between them, resulting in a weaker Coulomb repulsion, and, as a consequence, lower hydrophobic minimum. This effect can be harnessed to gain a better control over the overall shape of the potential. We limited ourselves to displacements that leave at least 3 Å between the charge and the inner boundary of the colloid so that situations where the charge is able to approach the wall at distances that permit chemical bonding are avoided. Figure 4b shows results obtained for colloid with R = 5 Å when using Cl− as counterions instead of F−. Both central and displaced charge configurations are shown. As shown in Figure 4a, the potential is purely repulsive for q = 3e and the central charge setup. But when the charge is displaced, a hydrophobic minimum is recovered. Shifting the charge, therefore, allows one to increase the maximum q, for which an SALR shape is observed. This results in a higher barrier for the same colloid size, as is clear from the comparison between the potentials for q = 2e, central charge, and q = 3e, shifted charge, models. As noted above, barriers can also be increased by increasing the size of the colloid. Figure 4c shows that the potential for R = 6 Å and q = 4e, which is qmax for this colloid, has a barrier of 5

Figure 5. Potential of mean force for the model with R = 7 Å and charge q = 3e displaced off center by 1 Å. Simulations were performed for Cl− counterions. The inset shows the transition temperature into the cluster crystal state as a function of the cluster size. For Nc ≥ 13, the crystal is predicted to be stable at room temperature.

R = 7 Å, Cl− are used as counterions and the charge of q = 3e is placed 1 Å away from the center. The potential has a minimum of −12 kJ/mol at rmin = 1.55 nm and a maximum of 5 kJ/mol at r rmax = 2.05 nm. The maximum to minimum ratio rmax = 1.32 is min

above the 1.29 threshold identified for the hard-sphere systems. An analytical extension was constructed for distances beyond rD = 2.2 nm, in which the screening constant κ = 1 nm−1 was used as an estimate suitable for the cluster crystal environment. The 2676

DOI: 10.1021/acs.jpcb.7b11662 J. Phys. Chem. B 2018, 122, 2669−2682

Article

The Journal of Physical Chemistry B

Figure 6. Colloid model with charge distributed over six particles. (a) Explanation of the geometry. Different colors denote different equivalent distances. The radius is 3.5 Å and the total charge q = 4e. Charges are elevated 1.3 Å above the colloid surface. (b) PMF for this model obtained in the presence of acetate salt at concentration corresponding to inverse screening length κ = 2.2 nm−1. The inset shows the temperature of the cluster fluid−cluster crystal transition. Crystals with cluster size 8 and higher are stable at temperatures below 290 K.

N*c . It also increases the maximum allowed charge qmax for which an SALR potential can be observed for the fixed size of the colloid. Alternatively, one can keep the charge fixed at q < qmax, but reduce the size of the colloid to achieve the same effect. This leads to a smaller total size of the system that needs to be simulated. The specifics of how exactly the charge is distributed define the overall shape of the potential. We note that the mutual repulsion between colloids in a close-contact configuration will go down when the charge is split into multiple particles. This property can be harnessed to adjust the depth of the minimum to a desired value in the same way as it is used for shifting of the charge to one side of the colloid. The screening and repulsion effects are subtle and interconnected, so finding a suitable model requires a thorough and careful search by simulations. 3.2.3.1. Model with Six Surface Charges. We studied a total of 68 systems with distributed charge. Five schemes of charge distribution were considered that differed in the number of subcharges, 2, 4, 6, 8, and 12, and the geometry of charge placement. Both symmetric and asymmetric charge configurations were examined. The radius of the colloid R was varied between 3 and 6 Å, whereas the total charge ranged from 3e to 6e. Two systems were seen to have a PMF that supports cluster crystal. The first has radius R = 5 Å and carries a total charge of 6e distributed over 12 particles that are placed on the surface on one side of the colloid. The second system is of smaller size, R = 3.5 Å, and carries a lower total charge of q = 4e. It has six sites with charge placed in a regular configuration on one side of the colloid. The geometrical details are explained in Figure 6a. Five charges are located at the corners of a pentagon, whereas the sixth charge is found at the center. All charges are at the same distance away from the center and elevated by 1.3 Å above the surface of the colloid. There are three distinct distances defining the geometry, as explained in Figure 6a. The results for the two systems are similar, so it suffices to discuss in detail just one of them. We will focus on the second system, which we will call 6P-system, as it is the smaller of the two and the easier one to simulate. The PMF of this system was extracted in a box with dimensions 5.5 × 3.1 × 3.1 nm3 filled with water and eight acetate ions. Analytical extension was built for r > rD = 1.5 nm and κ = 2.2 nm−1. The resulting potential, shown in Figure 6b,

same simulation protocol was applied as for the hard-sphere model system discussed in the preceding sections. First, simulated annealing was used to identify the ground-state configurations, which were seen to be the FCC cluster crystal. Second, replica-exchange simulations were performed for varying density to extract thermodynamic functions. The inset of Figure 5 shows the cluster fluid-to-cluster crystal transition temperature Tc as a function of the cluster size Nc. Except for a slightly higher cluster size and a plateau at large Nc, this dependence is identical to that seen for the model system in Figure 2f. Cluster crystals stable at room temperature are predicted for Nc ≥ Nc* = 13. According to the visual inspection, all colloids in such clusters can be considered located on the surface. Bulk colloids, by which we mean those that are surrounded by surface colloids from all directions, can be observed only for Nc ≥ 14. This agrees well with the observations made for atomic clusters.63,64 The simulation box for Nc = 13 contains 416 colloids, distributed among 32 equivalent clusters, and has the length of L = 20.13 nm. Unfortunately, explicit solvent simulations of a box of this size are prohibitive. As a consequence, more accurate tests of cluster crystals have to be carried out for a smaller system. 3.2.3. Model with Distributed Charge. There are two ways to reduce the size of the simulation box: (a) lower the number of colloids in the cluster and/or (b) consider smaller colloids. Below, we show that distributing the charge carried by the colloid across multiple sites helps to achieve both of these goals. Our tests for model potentials show that the size of the smallest cluster supported by the crystal is anticorrelated with the height of the maximum of the intercolloid potential, provided that the difference between the maximum and the minimum is kept ∼14 kJ/mol or more (data not shown). Delocalization of the charge is expected to reduce the amount of counterions bound to the colloid. For instance, for the model with R = 5 Å, q = 4e, F− counterions, and the central charge, one sees, on average, a charge of 1e bound to the colloid (situated 1.2 nm or less from the center) when simulated in a box with the side of 3.2 nm. The bound ionic charge drops to 0.44e when the colloid charge is broken into eight particles, each carrying 0.5e, and placed on the surface in a cubic configuration. Reduced screening increases the maximum in the effective potential between such particles, thus favoring smaller 2677

DOI: 10.1021/acs.jpcb.7b11662 J. Phys. Chem. B 2018, 122, 2669−2682

Article

The Journal of Physical Chemistry B

Figure 7. Results of the simulations in explicit solvent for 6P model. (a) Initial configuration in atomic detail. Subpanels demonstrate the distribution of acetate ions and water around colloids. (b) Q6 parameter computed in ∼10 ns trajectories for varying temperature. The time evolution for T = 330 K is illustrated in the four snapshots taken at notable points along the trajectory. The transition temperature can be estimated as Tc ≈ 320 K.

has a minimum of −10 kJ/mol at rmin = 0.87 nm and a maximum of 5.4 kJ/mol at rmax = 1.28 nm. The maximum to r minimum ratio rmax = 1.47 is well above the 1.29 threshold,

3.3. Cluster Lattice in Explicit Solvent. Explicit solvent simulations were conducted to test implicit-solvent predictions for the 6P model. To match the density observed in the effective-potential simulations, we use a system in which 32 clusters of size Nc = 9 are housed in a simulation box of length 9.97 nm. A box of slightly larger dimension, 10.22 nm, was considered initially, in which clusters were placed in an FCC configuration together with 1152 acetate ions at random positions. A constant-temperature simulation was conducted for 100 ps at T = 270 K with the positions of the colloids constrained to allow for the counterions to equilibrate. Another 100 ps run was conducted under constant-temperature, constant-pressure conditions with the constraints removed to allow the system achieve normal pressure of 1 atm. As a result of this latter equilibration run, the simulation box shrank to the desired size. Finally, a number of unconstrained constanttemperature simulations were performed for varying T to test the stability of the crystalline structure. The initial structure is presented in full atomic detail in Figure 7a. The colloid charges are found on the surface of the clusters and fully solvated by water. The counterions are distributed throughout the simulation box, but have a marked preference for the clusters. Water is distributed uniformly in the box. Figure 7b shows Q6 parameter computed for the locations of the cluster centroids at five temperatures, 270, 290, 310, 330, and 350 K. If we use 0.38 as the cutoff value to distinguish between crystalline and glassy states, then all trajectories fall into two categories: (a) those with T ≤ 310 K in which the crystal remains intact and (b) those with T ≥ 330 K in which the crystal is melting. A common feature in the trajectories that preserve the crystal is an initial transition period of ∼4 ns, over which the system reaches equilibrium. Q6 in that period drops from 0.56 to a

min

suggesting that the corresponding system should readily make cluster crystals. Indeed, simulations with the effective potential from Figure 6b find that the FCC configuration is the lowestenergy state. The transition temperature into this state Tc is shown in the inset of the figure as a function of the cluster size Nc. The crystal is stable at room temperature for Nc ≥ 8, where Tc ≥ 290. In addition to the FCC configurations, the simulations uncovered body-centered cubic (BCC) conformations, which have almost the same low potential energy. The prevalence of these configurations depends on the density and the temperature at which simulated annealing runs are conducted. The size of the simulation cell was also seen to have a non-negligible effect. Although it is difficult to determine exactly which configurations, FCC or BCC, represent the true free-energy minimum for the studied model, it is clear that polymorphic transitions at low temperatures cannot be excluded. A detailed analysis of the phase diagram will be the subject of a separate study. Here, we tested how the change of the ground-state configuration from FCC to BCC affects the results. Replicaexchange simulations were repeated for the new initial conformation, and new thermodynamic functions were obtained. The transition temperature was found to be the same as for FCC, or even slightly higher, depending on the density. We conclude, therefore, that the nature of the crystal phase in the cluster fluid-to-cluster crystal transition is not critical for the transition thermodynamics, at least for T ≈ 300 K. 2678

DOI: 10.1021/acs.jpcb.7b11662 J. Phys. Chem. B 2018, 122, 2669−2682

Article

The Journal of Physical Chemistry B

measured for a pair of colloids. At less than 10% error, the accuracy of this estimate is very high and can only be compared to the accuracy with which the temperature of the gas−liquid transition can be computed in the effective pair potential approximation, neglecting the nonadditivity effects, for the simplest of condensed matter systems, inert gases.73 In more complex systems, such as water, ionic liquids, or metals, the error is much larger and easily reaches 20% or more.74 More remarkably, the multibody effects are found to be very important even for the inert gases when it comes to the melting of solids instead of the gas−liquid transition. They contribute, for instance, to more than a 20% error in the melting temperature of argon, krypton, or xenon,75 or tip the balance toward the correct FCC structure in simulations of argon liquids.76 The liquid−crystal transition is also the subject of the present work, yet, making the results of our calculation all the more surprising, the accuracy of the transition temperature is much higher. Precisely, why the effective interactions turn out to be so accurate is not completely understood at the moment. Part of the reason is the design of the studied systems, in which the assembly of only small clusters is pursued. If all colloids in a cluster can be considered located on the surface, the effect of the multibody interactions is minimal. This is unlike the situation in the vapor−liquid transition, where particles in the dense liquid state face a radically different environment from what they experience in the dilute gas phase. However, in the liquid−solid transition, both phases are dense yet the multibody effects are still significant. Alternatively, the high accuracy of Tc could be due to the unorthodox shape of the interactions among colloids. Although individual colloids attract each other at short distances, clusters are separated by much longer distances and hence experience mostly repulsive interactions. These are the interactions that govern the transition temperature between the cluster fluid and cluster crystal phases. Is it possible that repulsive forces are not affected by the many-body effects to the same degree as are the attractive forces? In support of this thesis, binary mixtures of hard-sphere particles with different diameters, the ultimate models of repulsive interactions, are known to undergo a fluid− crystal transition that can be very accurately described by the effective pair interactions formulated for the larger particles alone.77 Whatever the real reason for the high accuracy of the effective interactions, they can be used as a guide in the search for actual physical systems experiencing a transition into the cluster crystal state. The features that are required for the crystal state are known and can be sought after in potentials produced experimentally, for instance, in small-angle scattering experiments with subsequent extraction of interparticle interactions, or in simulations, by computing PMFs for target compounds. There are good reasons to believe that such systems can be found. For the systems with charge concentrated in one particle, a good candidate are fullerenes.78 They come in about the right size range, R ∼ 10 Å, and can accommodate various atoms, ions, or clusters in their interior.79 Of particular interest could be metalofullerenes, the fullerenes with metal atom/ion embedded in their interior.80 There could be insurmountable challenges in designing such systems, not least because of the need to stabilize the charge and the location of the metals. The example of heme, which can cycle among several charge states

certain average value above 0.4, which depends on the temperature and then fluctuates around that average in the remaining ∼6 ns of the trajectory. The time scale of these fluctuations varies widely with temperature. At T = 270 K, the dynamics are very slow, but accelerate rapidly for T = 290 K and then again for T = 310 K. Judging from the figure, there is more than an order of magnitude difference in the relaxation times between the lowest and highest temperatures. Such a dramatic slowdown is puzzling and may suggest that the solvent in the system is close to a freezing point. We analyzed the behavior of water at T = 270 K and found no evidence of it being frozen or even supercooled. The mean-square displacement of water oxygen atoms recovers the expected diffusive regime within 2 ps for all temperatures. Furthermore, the extracted diffusion coefficient varies with temperature exactly in the same manner as for pure water so there is no anomaly at T = 270 K. The slow structural relaxation, therefore, must be attributed exclusively to the colloids themselves. One plausible explanation for slow dynamics is that free-energy barriers between different states that are all within the manifold of crystalline configurations, or the so-called basin of attraction on the free-energy landscape, grow rapidly with the decline of temperature. Such behavior is common among complex systems, in particular proteins. Among other things, it suggests a temperature regime in which the colloid subsystem is frozen, manifested in no colloid migration over available crystal sites, for instance, when the solvent remains in the liquid state. At T = 330 K, the relaxation time grows again. Snapshots at select moments in this trajectory are shown in the subpanels of Figure 7b. The order parameter reveals that the crystal melts within the first 1.5 ns. The snapshot at 2 ns, corresponding to Q6 = 0.35, clearly shows that the system has lost its crystalline structure. The loss is not irreversible, however. Later on at time 5.6 ns, the order parameter reaches a value of 0.43 and, as the corresponding snapshot illustrates, the order is partially restored. At later times, Q6 is always below 0.38 so the system remains disordered (see the snapshot at 11 ns). The long initial fluctuation, therefore, appears to be due to the system’s proximity to the phase transition. The relevant free-energy barrier here, unlike at the low temperatures, is the barrier separating the two stable phases, cluster crystal and cluster fluid. Because 310 K appears to be below the transition, and 330 K above it, we estimate that the transition temperature is Tc = 320 K. We note that it differs from the value of 290 K obtained for the same system interacting via effective potential by less than 10%. This is a remarkable accuracy given the severity of approximations involved in the implicit solvent treatment. At T = 350 K, the crystal melts within 1 ns of the simulation and never recovers. Simulations similar to those of the FCC lattice were conducted for the BCC lattice. The simulation box was 8.03 nm long and contained 16 clusters, each containing 9 colloids. In ∼10 ns runs, the crystal was stable at T = 270 and 290 K and melted at T = 350 K. The time traces of Q6 are very similar to those shown in Figure 7b and provide no new information. The crystal appears to be stable regardless of its specific structure.

4. DISCUSSION This work demonstrates that a combination of two interactions, (a) hydrophobic attraction and (b) electrostatic repulsion, is sufficient to induce a transition into the cluster crystal state. The transition temperature can be estimated at a low cost by implicit solvent model, in which effective interactions are 2679

DOI: 10.1021/acs.jpcb.7b11662 J. Phys. Chem. B 2018, 122, 2669−2682

Article

The Journal of Physical Chemistry B depending on its environment,81 shows just how precarious such a task could be. Fullerenes can also be functionalized with various charged groups.82 So, they also serve as a good candidate for the distributed-charge systems. Another example of such systems are proteins. We computed effective potentials for protein lysozyme in solutions with varying pH. At neutral pH, the interaction is attractive with a global minimum, precipitating liquid−liquid phase separation.83 Under acidic conditions, the protein picks up more charge and the global minimum is transformed into a local one, whereas the potential acquires a shape similar to that of type (a) potentials from Figure 2a. The minimum is at rmin = 26 Å and u(rmin) = 3.7 kJ/mol, whereas the maximum is at rmax = 36 Å and u(rmax) = 5.6 kJ/mol. Simulations show that this potential makes cluster crystals at low temperatures. For density ρ = 350 mg/mL, for instance, a lattice with four particles per site and 71 Å distance between nearest clusters is observed. This density is easily achieved experimentally; at pH 6, it leads to a liquid−liquid phase separation.83,84 The transition temperature Tc into the crystalline state, however, is only about 70 K. With increasing density, columnar structures begin to emerge instead of the cluster crystal. So, this strategy cannot be used to increase Tc. As Figure 2a suggests, for transitions to occur at room temperature, the barrier in the potential has to go up roughly three times. Whether this can be achieved for lysozyme by further lowering pH is not clear. Perhaps, mutations increasing the number of charged groups on the protein surface would help. In any case, the example of lysozyme shows that proteins make good candidates for cluster-crystal-making systems. Inorganic macroions, such as polyoxometalates (POMs)85 and metal chalcogenide clusters,86 represent another class of molecules that can hypothetically form cluster crystals. A major advantage of these compounds over proteins is the presence of a hard core with defined geometry, which can be easily controlled through synthetic approaches with atomic precision. However, POMs have a large amount of surface oxygen atoms, which can form specific hydrogen bonds with water molecules87 so that they should be protected with organic coating to maintain their hydrophobic nature. Monodisperse nanocrystals (NCs) present an extension of inorganic nanoclusters to larger sizes. Their hydrophobic inorganic core can be covered with various ligands, allowing precise control of their charge and hydrophobic properties. The recent investigations of Kotov et al.88 show evidence supporting the formation of self-limiting colloidal clusters from different NCs capped with negatively charged citrate or poly(styrene sulfonate) polymer. An unexpected result is that monodisperse clusters can be obtained from originally polydisperse nanoparticles. In addition, a micrometer-sized colloidal crystal was prepared from these “supraparticles”. However, it is unlikely that these clusters and corresponding cluster lattices are in equilibrium because both citrate and poly(styrene sulfonate) can act as bridging ligands between neighboring NCs forming strong bonds with them.

■

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Computational resources of the research cluster managed by the Institute for Condensed Matter Physics, National Academy of Sciences of Ukraine, and the Ukrainian National Grid were used for the completion of this work. W.C. was supported by US National Science Foundation (grant DMS-1764187).

■

REFERENCES

(1) Li, M.; Ober, C. K. Block Copolymer Patterns and Templates. Mater. Today 2006, 9, 30−39. (2) Klein, W.; Gould, H.; Ramos, R. A.; Clejan, I.; Melcuk, A. I. Repulsive Potentials, Clumps and the Metastable Glass Phase. Phys. A 1994, 205, 738−746. (3) Lenz, D. A.; Mladek, B. M.; Likos, C. N.; Blaak, R. Thermodynamic Stability and Structural Properties of Cluster Crystals Formed by Amphiphilic Dendrimers. J. Chem. Phys. 2016, 144, No. 204901. (4) Likos, C. N.; Lang, A.; Watzlawek, M.; Löwen, H. Criterion for Determining Clustering Versus Reentrant Melting Behavior for Bounded Interaction Potentials. Phys. Rev. E 2001, 63, No. 031206. (5) Wilding, N. B.; Sollich, P. Demixing Cascades in Cluster Crystals. J. Chem. Phys. 2014, 141, No. 094903. (6) Likos, C. N. Colloidal Interactions: From Effective Potentials to Structure. Riv. Nuovo Cimento 2014, 37, 125−180. (7) Pàmies, J. C.; Cacciuto, A.; Frenkel, D. Phase Diagram of Hertzian Spheres. J. Chem. Phys. 2009, 131, No. 044514. (8) Mladek, B. M.; Charbonneau, P.; Frenkel, D. Phase Coexistence of Cluster Crystals: Beyond the Gibbs Phase Rule. Phys. Rev. Lett. 2007, 99, No. 235702. (9) Neuhaus, T.; Likos, C. N. Phonon Dispersions of Cluster Crystals. J. Phys.: Condens. Matter 2011, 23, No. 234112. (10) Moreno, A. J.; Likos, C. N. Diffusion and Relaxation Dynamics in Cluster Crystals. Phys. Rev. Lett. 2007, 99, No. 107801. (11) Nikoubashman, A.; Kahl, G.; Likos, C. N. Cluster Crystals under Shear. Phys. Rev. Lett. 2011, 107, No. 068302. (12) Slimani, M. Z.; Bacova, P.; Bernabei, M.; Narros, A.; Likos, C. N.; Moreno, A. J. Cluster Glasses of Semiflexible Ring Polymers. ACS Macro Lett. 2014, 3, 611−616. (13) Sciortino, F.; Zaccarelli, E. Computational Materials Science Soft Heaps and Clumpy Crystals. Nature 2013, 493, 30−31. (14) Lenz, D. A.; Blaak, R.; Likos, C. N.; Mladek, B. M. Microscopically Resolved Simulations Prove the Existence of Soft Cluster Crystals. Phys. Rev. Lett. 2012, 109, No. 228301. (15) Ziherl, P.; Kamien, R. D. From Lumps to Lattices: Crystallized Clusters Made Simple. J. Phys. Chem. B 2011, 115, 7200−7205. (16) Košmrlj, A.; Pauschenwein, G. J.; Kahl, G.; Ziherl, P. Continuum Theory for Cluster Morphologies of Soft Colloids. J. Phys. Chem. B 2011, 115, 7206−7217. (17) Salgado-Blanco, D.; Mendoza, C. I. Self-Assembly of Anisotropic Soft Particles in Two Dimensions. Eur. Phys. J. E 2013, 36, 38. (18) Zhuang, Y.; Charbonneau, P. Recent Advances in the Theory and Simulation of Model Colloidal Microphase Formers. J. Phys. Chem. B 2016, 120, 7775−7782. (19) Noro, M. G.; Frenkel, D. Extended Corresponding-States Behavior for Particles with Variable Range Attractions. J. Chem. Phys. 2000, 113, 2941−2944. (20) Kendrick, G. F.; Sluckin, T. J.; Grimson, M. J. Macrocrystal Phases in Charged Colloidal Suspensions. Europhys. Lett. 1988, 6, 567−572. (21) Seul, M.; Andelman, D. Domain Shapes and Patterns: the Phenomenology of Modulated Phases. Science 1995, 267, 476−483. (22) Sear, R. P.; Gelbart, W. M. Microphase Separation Versus the Vapor−Liquid Transition in Systems of Spherical Particles. J. Chem. Phys. 1999, 110, 4582−4588.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

A. Baumketner: 0000-0003-2726-931X 2680

DOI: 10.1021/acs.jpcb.7b11662 J. Phys. Chem. B 2018, 122, 2669−2682

Article

The Journal of Physical Chemistry B

(46) Kaminski, G. A.; Friesner, R. A.; Tirado-Rives, J.; Jorgensen, W. L. Evaluation and Reparametrization of the OPLS-AA Force Field for Proteins Via Comparison with Accurate Quantum Chemical Calculations on Peptides. J. Phys. Chem. B 2001, 105, 6474−6487. (47) Kaminski, G.; Duffy, E. M.; Matsui, T.; Jorgensen, W. L. FreeEnergies of Hydration and Pure Liquid Properties of Hydrocarbons from the OPLS All-Atom Model. J. Phys. Chem. 1994, 98, 13077− 13082. (48) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Single Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (49) Torrie, G. M.; Valleau, J. P. Non-Physical Sampling Distributions in Monte-Carlo Free-Energy Estimation: Umbrella Sampling. J. Comput. Phys. 1977, 23, 187−199. (50) Ferrenberg, A. M.; Swendsen, R. H. New Monte-Carlo Technique for Studying Phase-Transitions. Phys. Rev. Lett. 1988, 61, No. 2635. (51) Chandler, D. Interfaces and the Driving Force of Hydrophobic Assembly. Nature 2005, 437, 640−647. (52) Sharp, K. A.; Nicholls, A.; Fine, R. F.; Honig, B. Reconciling the Magnitude of the Microscopic and Macroscopic Hydrophobic Effects. Science 1991, 252, 106−109. (53) Lee, B.; Richards, F. M. Interpretation of Protein Structures: Estimation of Static Accessibility. J. Mol. Biol. 1971, 55, 379−389. (54) Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 2005, 26, 1701−1718. (55) Chandrasekhar, J.; Spellmeyer, D. C.; Jorgensen, W. L. Energy Component Analysis for Dilute Aqueous Solutions of Li+, Na+, F−, and Cl− Ions. J. Am. Chem. Soc. 1984, 106, 903−910. (56) McDonald, N. A.; Duffy, E. M.; Jorgensen, W. L. Monte Carlo Investigations of Selective Anion Complexation by a Bis(Phenylurea)P-Tert-Butylcalix[4]Arene. J. Am. Chem. Soc. 1998, 120, 5104−5111. (57) Ȧ qvist, J. Ion Water Interaction Potentials Derived from FreeEnergy Perturbation Simulations. J. Phys. Chem. 1990, 94, 8021−8024. (58) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (59) Miyamoto, S.; Kollman, P. A. Settle: an Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models. J. Comput. Chem. 1992, 13, 952−962. (60) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463−1472. (61) Nosé, S. Constant Temperature Molecular-Dynamics Methods. Prog. Theor. Phys. Suppl. 1991, 1−46. (62) Roux, B. Implicit Solvent Models. In Computational Biochemistry and Biophysics; CRC Press, 2001. (63) Sun, Y.; Zhang, M.; Fournier, R. Periodic Trends in the Geometric Structures of 13-Atom Metal Clusters. Phys. Rev. B 2008, 77, No. 075435. (64) Fernando, A.; Weerawardene, K. L. D. M.; Karimova, N. V.; Aikens, C. M. Quantum Mechanical Studies of Large Metal, Metal Oxide, and Metal Chalcogenide Nanoparticles and Clusters. Chem. Rev. 2015, 115, 6112−6216. (65) Sciortino, F.; Tartaglia, P.; Zaccarelli, E. One-Dimensional Cluster Growth and Branching Gels in Colloidal Systems with ShortRange Depletion Attraction and Screened Electrostatic Repulsion. J. Phys. Chem. B 2005, 109, 21942−21953. (66) Steinhardt, P. J.; Nelson, D. R.; Ronchetti, M. BondOrientational Order in Liquids and Glasses. Phys. Rev. B 1983, 28, No. 784. (67) Leocmach, M.; Tanaka, H. Roles of Icosahedral and Crystal-Like Order in the Hard Spheres Glass Transition. Nat. Commun. 2012, 3, No. 974. (68) Götzelmann, B.; Evans, R.; Dietrich, S. Depletion Forces in Fluids. Phys. Rev. E 1998, 57, No. 6785.

(23) Baumketner, A.; Cai, W. Equilibrium Clusters in Suspensions of Colloids Interacting Via Potentials with a Local Minimum. Condens. Matter Phys. 2016, 19, No. 13605. (24) Ciach, A. Universal Sequence of Ordered Structures Obtained from Mesoscopic Description of Self-Assembly. Phys. Rev. E 2008, 78, No. 061505. (25) Sciortino, F.; Mossa, S.; Zaccarelli, E.; Tartaglia, P. Equilibrium Cluster Phases and Low-Density Arrested Disordered States: The Role of Short-Range Attraction and Long-Range Repulsion. Phys. Rev. Lett. 2004, 93, No. 055701. (26) Archer, A. J.; Pini, D.; Evans, R.; Reatto, L. Model Colloidal Fluid with Competing Interactions: Bulk and Interfacial Properties. J. Chem. Phys. 2007, 126, No. 014104. (27) Valadez-Pérez, N. E.; Castañeda-Priego, R.; Liu, Y. Percolation in Colloidal Systems with Competing Interactions: The Role of LongRange Repulsion. RSC Adv. 2013, 3, 25110−25119. (28) Godfrin, P. D.; Castañeda-Priego, R.; Liu, Y.; Wagner, N. J. Intermediate Range Order and Structure in Colloidal Dispersions with Competing Interactions. J. Chem. Phys. 2013, 139, No. 154904. (29) Cao, X. J.; Cummins, H. Z.; Morris, J. F. Hydrodynamic and Interparticle Potential Effects on Aggregation of Colloidal Particles. J. Colloid Interface Sci. 2012, 368, 86−96. (30) Bomont, J.-M.; Bretonnet, J.-L.; Costa, D. Temperature Study of Cluster Formation in Two-Yukawa Fluids. J. Chem. Phys. 2010, 132, No. 184508. (31) Richardi, J. One-Dimensional Assemblies of Charged Nanoparticles in Water: A Simulation Study. J. Chem. Phys. 2009, 130, No. 044701. (32) Archer, A. J.; Wilding, N. B. Phase Behavior of a Fluid with Competing Attractive and Repulsive Interactions. Phys. Rev. E 2007, 76, No. 031501. (33) Toledano, J. C. F.; Sciortino, F.; Zaccarelli, E. Colloidal Systems with Competing Interactions: From an Arrested Repulsive Cluster Phase to a Gel. Soft Matter 2009, 5, 2390−2398. (34) Liu, Y.; Porcar, L.; Chen, J.; Chen, W.-R.; Falus, P.; Faraone, A.; Fratini, E.; Hong, K.; Baglioni, P. Lysozyme Protein Solution with an Intermediate Range Order Structure. J. Phys. Chem. B 2011, 115, 7238−7247. (35) de Candia, A.; Del Gado, E.; Fierro, A.; Sator, N.; Tarzia, M.; Coniglio, A. Columnar and Lamellar Phases in Attractive Colloidal Systems. Phys. Rev. E 2006, 74, No. 010403. (36) Coniglio, A.; de Candia, A.; Fierro, A. Modulated Phases and Structural Arrest in Colloidal Systems with Competing Interactions. Mol. Phys. 2011, 109, 2981−2987. (37) Zhuang, Y.; Charbonneau, P. Equilibrium Phase Behavior of the Square-Well Linear Microphase-Forming Model. J. Phys. Chem. B 2016, 120, 6178−6188. (38) Zhuang, Y.; Zhang, K.; Charbonneau, P. Equilibrium Phase Behavior of a Continuous-Space Microphase Former. Phys. Rev. Lett. 2016, 116, No. 098301. (39) Almarza, N. G.; Pȩkalski, J.; Ciach, A. Effects of Confinement on Pattern Formation in Two Dimensional Systems with Competing Interactions. Soft Matter 2016, 12, 7551−7563. (40) Swendsen, R. H.; Wang, J. S. Replica Monte-Carlo Simulation of Spin-Glasses. Phys. Rev. Lett. 1986, 57, No. 2607. (41) Hukushima, K.; Nemoto, K. Exchange Monte Carlo Method and Application to Spin Glass Simulations. J. Phys. Soc. Jpn. 1996, 65, 1604−1608. (42) Hansmann, U. H. E. Parallel Tempering Algorithm for Conformational Studies of Biological Molecules. Chem. Phys. Lett. 1997, 281, 140−150. (43) Sugita, Y.; Okamoto, Y. Replica-Exchange Molecular Dynamics Method for Protein Folding. Chem. Phys. Lett. 1999, 314, 141−151. (44) Ferrenberg, A. M.; Swendsen, R. H. Optimized Monte-Carlo Data-Analysis. Phys. Rev. Lett. 1989, 63, No. 1195. (45) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. The Weighted Histogram Analysis Method for FreeEnergy Calculations on Biomolecules. 1. The Method. J. Comput. Chem. 1992, 13, 1011−1021. 2681

DOI: 10.1021/acs.jpcb.7b11662 J. Phys. Chem. B 2018, 122, 2669−2682

Article

The Journal of Physical Chemistry B (69) Chung, S.-W.; Markovich, G.; Heath, J. R. Fabrication and Alignment of Wires in Two Dimensions. J. Phys. Chem. B 1998, 102, 6685−6687. (70) Sear, R. P.; Chung, S.-W.; Markovich, G.; Gelbart, W. M.; Heath, J. R. Spontaneous Patterning of Quantum Dots at the Air− Water Interface. Phys. Rev. E 1999, 59, No. R6255. (71) Alberts, B.; Johnson, A.; Lewis, J.; Morgan, D.; Raff, M.; Roberts, K.; Walter, P. Molecular Biology of the Cell; Garland Science, 2014. (72) Fowler, R. H.; Guggenheim, E. A. Statistical Thermodynamics: A Version of Statistical Mechanics for Students of Physics and Chemistry; University Press: Cambridge, England, 1965. (73) Bukowski, R.; Szalewicz, K. Complete Ab Initio Three-Body Nonadditive Potential in Monte Carlo Simulations of Vapor−Liquid Equilibria and Pure Phases of Argon. J. Chem. Phys. 2001, 114, 9518− 9531. (74) Piela, L. Intermolecular Interactions A2. In Ideas of Quantum Chemistry; Elsevier: Amsterdam, 2007; Chapter 13, pp 681−761. (75) Wang, L.; Sadus, R. J. Three-Body Interactions and Solid− Liquid Phase Equilibria: Application of a Molecular Dynamics Algorithm. Phys. Rev. E 2006, 74, No. 031203. (76) Lotrich, V. F.; Szalewicz, K. Three-Body Contribution to Binding Energy of Solid Argon and Analysis of Crystal Structure. Phys. Rev. Lett. 1997, 79, No. 1301. (77) Dijkstra, M.; van Roij, R.; Evans, R. Direct Simulation of the Phase Behavior of Binary Hard-Sphere Mixtures: Test of the Depletion Potential Description. Phys. Rev. Lett. 1999, 82, No. 117. (78) Kroto, H. W.; Heath, J. R.; O’Brien, S. C.; Curl, R. F.; Smalley, R. E. C-60: Buckminsterfullerene. Nature 1985, 318, 162−163. (79) del Carmen Gimenez-Lopez, M.; Chuvilin, A.; Kaiser, U.; Khlobystov, A. N. Functionalised Endohedral Fullerenes in SingleWalled Carbon Nanotubes. Chem. Commun. 2011, 47, 2116−2118. (80) Chai, Y.; Guo, T.; Jin, C. M.; Haufler, R. E.; Chibante, L. P. F.; Fure, J.; Wang, L. H.; Alford, J. M.; Smalley, R. E. Fullerenes with Metals Inside. J. Phys. Chem. 1991, 95, 7564−7568. (81) Li, Z. The Secret Life of Heme in Regulating Diverse Biological Processes. In Heme Biology; World Scientific: Amsterdam, 2011. (82) Periya, V. K.; Koike, I.; Kitamura, Y.; Iwamatsu, S.-i.; Murata, S. Hydrophilic [60]Fullerene Carboxylic Acid Derivatives Retaining the Original 60π Electronic System. Tetrahedron Lett. 2004, 45, 8311− 8313. (83) Baumketner, A.; Melnyk, R.; Holovko, M. F.; Cai, W.; Costa, D.; Caccamo, C. Softness and Non-Spherical Shape Define the Phase Behavior and the Structural Properties of Lysozyme in Aqueous Solutions. J. Chem. Phys. 2016, 144, No. 015103. (84) Kastelic, M.; Kalyuzhnyi, Y. V.; Hribar-Lee, B.; Dill, K. A.; Vlachy, V. Protein Aggregation in Salt Solutions. Proc. Natl. Acad. Sci. U.S.A. 2015, 112, 6766−6770. (85) Yin, P.; Li, D.; Liu, T. Solution Behaviors and Self-Assembly of Polyoxometalates as Models of Macroions and Amphiphilic Polyoxometalate−Organic Hybrids as Novel Surfactants. Chem. Soc. Rev. 2012, 41, 7368−7383. (86) Dance, I.; Fisher, K.; Wachter, J.; Fenske, D.; Ohmer, J.; Hachgenei, J.; Merzweiler, K.; Muller, A.; Jostes, R.; Karlin, F. C. K. Progress in Inorganic Chemistry; John Wiley & Sons: New York, 1994; Vol. 41, p 637. (87) Chaumont, A.; Wipff, G. Do Keggin Anions Repulse Each Other in Solution? The Effect of Solvent, Counterions and Ion Representation Investigated by Free Energy (PMF) Simulations. C. R. Chim. 2012, 15, 107−117. (88) Xia, Y.; Nguyen, T. D.; Yang, M.; Lee, B.; Santos, A.; Podsiadlo, P.; Tang, Z.; Glotzer, S. C.; Kotov, N. A. Self-Assembly of SelfLimiting Monodisperse Supraparticles from Polydisperse Nanoparticles. Nat. Nanotechnol. 2011, 6, 580−587.

2682

DOI: 10.1021/acs.jpcb.7b11662 J. Phys. Chem. B 2018, 122, 2669−2682