Low-Temperature Crystal Structures of the Hard Core Square Shoulder Model

In many cases, the stability of complex structures in colloidal systems is enhanced by a competition between different length scales. Inspired by recent experiments on nanoparticles coated with polymers, we use Monte Carlo simulations to explore the types of crystal structures that can form in a simple hard-core square shoulder model that explicitly incorporates two favored distances between the particles. To this end, we combine Monte Carlo-based crystal structure finding algorithms with free energies obtained using a mean-field cell theory approach, and draw phase diagrams for two different values of the square shoulder width as a function of the density and temperature. Moreover, we map out the zero-temperature phase diagram for a broad range of shoulder widths. Our results show the stability of a rich variety of crystal phases, such as body-centered orthogonal (BCO) lattices not previously considered for the square shoulder model.

Here we introduce a variation on the "floppy box" Monte Carlo (FBMC) crystal prediction approach proposed in Ref. [1,2]. In the FBMC approach, candidate crystal structures are predicted by performing a Monte Carlo simulation for a small number of particles (typically N 20). The particles are contained in a periodic simulation box spanned by three vectors that can fluctuate in both length and direction. Due to the small number of particles, the system can sample a large number of crystal phases with differently shaped unit cells in a short time, allowing it to rapidly find entropically and energetically favorable structures. These structures can then be tested for stability using simulations on a larger scale or theoretical approximations.
The efficiency of simulations of small boxes with a variable shape is highly dependent on the degree of distortion of the simulation box. In particular, if the box is flattened in any direction, particles can interact with periodic images of themselves and other particles that are several boxes away. To keep the box distortion low, the floppy box approach typically makes use of a lattice reduction technique [3], which redescribes a given unit cell in terms of a linear combination of its box vectors.
In order to simplify this extra step in the simulation, we propose a variation on the FBMC approach which makes use of shifted boundary conditions. In this shiftedboundary Monte Carlo (SBMC) approach the unit cell is defined by three perpendicular box vectors where L x,y,z are the box lengths along the x, y, and zaxis. Periodic boundary conditions are then shifted such that the origin of the next periodic image of the box in the three directions are given by: where δ ij is the relative shift along the i-axis associated with moving one periodic image in the j-direction. Note that the vectors p x,y,z span any unit cell shaped like a parallelepiped. In our simulations, we now simply vary the values of L x,y,z and δ ij independently in our Monte Carlo simulation, sampling different crystal structures in a way that is equivalent to the FBMC method. Effectively, the only change we have made thusfar is the requirement that one of the vectors spanning the periodic unit cell lies along the x-axis, and one lies in the xy-plane, eliminating the (irrelevant) rotational freedom of the full unit cell.
However, in this representation, lattice reduction can be done by a simple redefinition of the shifting parameters δ ij such that they each lie between −0.5 and 0.5. In particular, starting from set of δ ij , the reduced shifting parameters are given by where [·] denotes rounding to the nearest integer, and is a modulo function which always returns a number between −0.5 and 0.5. It is important to note that the reduction of the shifting parameters here does not completely prevent the formation of flat boxes, analogous to the lattice reduction in the FBMC approach. For example, particles can still line up into a widely spaced column in which the particles interact with periodic images several boxes away. Following Ref. [1], we prevent this by implementing a minimum axis length for all three directions, requiring that L x,y,z > L min = 0.75σ, with σ the size of the hard core of our particles. No angular restrictions are required.
We store our particle coordinates as scaled coordinates s, such that they can be converted to real-space coordinates via As a result, when performing a Monte Carlo move which changes one of the parameters L x,y,z , particle positions scale along with the volume change, and we apply the standard Monte Carlo acceptance rule for volume rescaling (see e.g. Ref. [2]). When performing a move which changes one of the shifting parameters δ ij , the volume remains constant, and acceptance is purely based on the change in energy of the system.
To determine which periodic image particles we should take into account when calculating interactions, we make use of our minimum axis length L min . Under the condition that −0.5 ≤ δ ij ≤ 0.5, the origin of an image of the simulation box that is n boxes away along the x, y, or z direction is at least a distance nL min away from the origin of the simulation box. As a result, any particle in this image box is at least a distance (n − 1)L min . Hence, for particles with a maximum interaction range σ + δ, we at most have to consider image particles up to (σ + δ)/L min + 1 boxes away in each direction. We loop over all of these boxes, and determine the actual minimum distance between a point in the image box and a point in the central box. We then take into account image particles in each image box where this distance is less than the interaction range.
We have confirmed that the FBMC and SBMC approaches largely predict the same structures for our square-shoulder systems with interaction ranges δ/σ = 0.15 and 0.20. Although the results are not expected to be identical, both methods find the same set of stable zero-temperature structures, and show strong similarity in the frequencies at which different structures are observed.
It is important to note that the modification to the FBMC approach proposed here is neither significantly faster or more effective for the purpose of predicting candidate crystal structures. However, implementation of the SBMC approach avoids the lattice reduction technique and simplifying the constraints for avoiding extreme distortion of the simulation box. As a result, the SBMC approach is considerably easier to implement.

CELL THEORY USING VORONOI CELLS FOR INTERACTING PARTICLES
Cell theory [4,5] provides a simple route for calculating a mean-field approximation for the free energy of a crystalline solid. In cell theory, the partition function of a single particle in the crystal is approximated by assuming that all other particles are located at their lattice site. The partition function of the particle under consideration is then written as: where β = 1/k B T with k B Boltzmann's constant and T the temperature, and u(r) is the energy of the particle at position r. The thermal wavelength Λ has no effect on the phase behavior, and can therefore freely be chosen to equal the particle diameter σ. In the case of hard spherical particles, this partition function reduces to the free volume available to the central particle (divided by Λ 3 ). This is commonly approximated by assuming that this volume is a polyhedron, formed by taking the Wigner-Seitz cell of the particle, and moving each of its faces inward to a distance r i − σ from the lattice position of the central particle, where r i is the distance to neighbor i in the perfect lattice [6]. Although this slightly underestimates the free volume available to the central particle, the effects of this approximation are on the same order as the approximations made in cell theory, and the volume of the resulting polyhedron is significantly easier to calculate than the "true" free volume.
Here, we extend this approach to particles interacting with square shoulder repulsions, with hard core diameter σ, shoulder width δ, and shoulder height . In this case, we consider for each of the N neighbors of the central particle two possible planes delimiting the free volume of the central particle: one at distance σ from the neighbor position (representing its hard-core repulsion), and one at distance σ + δ from the neighbor position (representing its shoulder), as illustrated in Fig. 1a. Clearly, all of the outer planes again form a polyhedron which approximates the free volume of the central particles if the square-shoulder repulsion was absent (blue in Fig. 1b). Similarly, the innermost polyhedron represents the volume V ∅ available to the particle where it interacts with no other particles, i.e. its energy u = 0. In total, we can construct 2 N polyhedra, each associated with a different combination of particle interactions. For example, the red polyhedron drawn in Fig. 1c shows the volume V {5} available to the particle if it at most interacts with neighbor 5. We can calculate the subvolume associated with the case where the particle interacts only with neighbor 5 (shaded blue in Fig. 1c) Using the same approach, we can calculate e.g. the subvolume V {5,6} (Fig. 1d) where the particle interacts with both particle 5 and 6 via This approach can be straightforwardly repeated until the subvolume for each set of interactions is known. The total single-particle partition function can then be written as: where δV Si denotes the subvolume for interacting with particles in the set S i of neighbors, and u Si is the energy of that state, given by times the length of S i . Additionally, u ref is the energy of the particle in the central position.
To calculate these volumes, we construct 2 N polyhedra and calculate their volume V i using the Voro++ library [7]. Any cell that does not contain the central lattice position is assumed to have volume V i = 0. This only occurs if the ground state energy of the crystal cell under consideration has an energy u ref > 0. Effectively, this means that we ignore the energy change associated with moving away from a neighbor the particle interacts with at its ideal lattice site. This additional approximation is justified by the observation that if these positions contribute meaningfully to the free energy, the cell contains local energy minima which do not correspond to the central lattice position, implying that the particle has multiple favored positions within its cell. Since this conflicts with the mean-field assumption that all particles are (on average) at their lattice site, cell theory is already likely to break down in these cases. Moreover, a crystal structure where particles can lower their energy by moving away from their lattice site is not likely to be stable. After calculating all volumes, we calculate the associated subvolumes by subtracting from each volume V i all subvolumes that correspond to neighbor sets that are strict subsets of neighbor set i. If we set δV 0 = V 0 , we can write this as: To test the effect of the additional approximations made in this method, we compared the results of this approach to cell theory free energies calculated using Monte Carlo integration of the single-particle partition function in Eq. 5, performed by repeated insertion of the central particle into a small insertion volume around its central position. The results are in good agreement within the deviations normally expected from cell theory. As an example, we plot in Fig. 2 a comparison of the free energy of a body-centered cubic crystal for interaction range δ/σ = 0.15 as a function of density at fixed temperature, calculated using both random insertion and the Voronoi cell approach.