From square-well to Janus: Improved algorithm for integral equation theory and comparison with thermodynamic perturbation theory within the Kern-Frenkel model

From square-well to Janus: Improved algorithm for integral equation theory and comparison with thermodynamic perturbation theory within the Kern-Frenkel model

Achille Giacometti Dipartimento di Scienze Molecolari e Nanosistemi, Università Ca’ Foscari Venezia, Calle Larga S. Marta DD2137, I-30123 Venezia, Italy    Christoph Gögelein Max-Planck-Institute for Dynamics and Self-Organization, Göttingen, Germany    Fred Lado Department of Physics, North Carolina State University, Raleigh, North Carolina 27695-8202    Francesco Sciortino Dipartimento di Fisica and CNR-SOFT, Università di Roma La Sapienza, Piazzale A. Moro 2, 00185 Roma, Italy    Silvano Ferrari Institut für Theoretische Physik and Center for Computational Materials Science, Technische Universität Wien, Wiedner Hauptstraße 8-10/136, A-1040 Wien, Austria    Giorgio Pastore Dipartimento di Fisica dell’ Università di Trieste and CNR-IOM, Strada Costiera 11, 34151 Trieste, Italy
July 7, 2019

Building upon past work on the phase diagram of Janus fluids [Sciortino et al., Phys. Rev. Lett. 103, 237801 (2009)], we perform a detailed study of integral equation theory of the Kern-Frenkel potential with coverage that is tuned from the isotropic square-well fluid to the Janus limit. An improved algorithm for the reference hypernetted-chain (RHNC) equation for this problem is implemented that significantly extends the range of applicability of RHNC. Results for both structure and thermodynamics are presented and compared with numerical simulations. Unlike previous attempts, this algorithm is shown to be stable down to the Janus limit, thus paving the way for analyzing the frustration mechanism characteristic of the gas-liquid transition in the Janus system. The results are also compared with Barker-Henderson thermodynamic perturbation theory on the same model. We then discuss the pros and cons of both approaches within a unified treatment. On balance, RHNC integral equation theory, even with an isotropic hard-sphere reference system, is found to be a good compromise between accuracy of the results, computational effort, and uniform quality to tackle self-assembly processes in patchy colloids of complex nature. Further improvement in RHNC however clearly requires an anisotropic reference bridge function.

Integral equation theory, Janus fluid, phase diagrams

I Introduction

Stimulated by recent advances in chemical syntheses of colloidal particles with different forms and functionalities, Walther09 (); Pawar10 () theoretical approaches have made significant progress in the last few years. Patchy colloids Glotzer04 (); Glotzer07 () in particular, having their surfaces decorated with different functionalities (e.g., solvophobic in opposition to solvophilic moieties), appear to combine the possibility of obtaining a large number of targeted structures, on the one hand, along with the possibility of local rearrangements, on the other hand, that represent the optimal trade-off for engineering self-assembly processes at mesoscopic scales. Whitesides02 ()

While direct comparison of theory with experiment still relies heavily on extensive numerical simulations that constitute today the main theoretical tool, given their virtually exact predictions, the heavy computational effort imposed by the anisotropic nature of patchy interactions (see e.g. Refs. Sciortino09, ; Sciortino10, ) has stimulated attempts to find approximate, yet reliable, alternative methods that can provide semi-quantitative estimates within a modest amount of computer time.

Two of these methods with established roles in liquid state studies Gray84 (); Hansen86 () are integral equation theory and thermodynamic perturbation theory. The main aim of integral equation theory is the computation of the pair correlation function, from which one can derive all thermodynamic and structural quantities. In order to perform practical computations, one is forced to introduce here an approximation into the exact relation between pair potential and pair distribution function, i.e. selecting a closure equation. In thermodynamic perturbation theory, on the other hand, the free energy of the system can be computed as a perturbation series of terms, provided the free energy and many-particle distribution functions of a reference systems are known. Usually the expansion is approximated by the truncation of the infinite series to the few terms that can be evaluated.

In the present paper, we discuss the performances of both methods when applied to a particular model, the Kern-Frenkel potential Kern03 (); Chapman88 () for patchy colloids, that has recently proven very useful within this anisotropic framework. Building upon previous work, Giacometti09a (); Giacometti10 (); Giacometti09b (); Gogelein12 () we compare the performance of a specific integral equation closure, the reference hypernetted-chain (RHNC), Lado73 (); Rosenfeld79 () and of a specific thermodynamic perturbation theory, devised by Barker and Henderson (TPT-BH), Barker67 (); Barker76 () on the single-patch Kern-Frenkel potential. In the case of the RHNC integral equation, generalized for molecular fluids, Lado82a (); Lado82b () we additionally present an improved algorithm allowing us to reach the limit of equal solvophobic-solvophilic composition, known as the Janus limit, that was not reachable with the original algorithm presented in Ref. Giacometti09a, .

The remainder of the paper is organized as follows. In Sec. II, we briefly recall the Kern-Frenkel model, while in Sec. III and Sec. V we review the application to this problem of the RHNC integral equation approach of Ref. Giacometti09a, and the TPT-BH of Ref. Gogelein12, . The improved algorithm for RHNC is described in Sec. IV and a detailed comparison of the performance of the two methods in contrast to numerical simulations is provided in Sec. VI. Section VII completes the paper with some conclusions and perspectives.

Ii The Kern-Frenkel model

The model for patchy interactions in colloids that we study here is due to Kern and Frenkel, Kern03 () an elaboration of the original model by Chapman et al. Chapman88 () They consider a fluid of hard spheres where the surface of each sphere is divided into two parts having square-well and hard-sphere character, the first mimicking a solvophobic region, the second a solvophilic region, within an implicit solvent description. Because of the azimuthal symmetry, the angular width of the solvophobic region is described by a single polar angle that becomes equal to in the even-division case (the Janus limit).

The positions of the particles in volume are given by a set of vectors , with , while the angular orientation of each square-well patch on a sphere surface is identified by unit vector . Finally, the direction connecting the centers of spheres and is characterized by unit vector . Figure 1 depicts the situation in the case of the Janus limit.

Figure 1: The one-patch Kern-Frenkel model, where is the direction joining the two centers and the orientations of the patches are specified by unit vectors and . The present configuration depicts the Janus limit.

Thus, two spheres of diameter attract each other via a square-well potential of width and depth , if the directions of the patch on each sphere are within a solid angle defined by and their relative distance lies within the range of the attractive well, and repel each other as hard spheres otherwise. As the system is still translationally invariant, the pair potential depends upon the difference , rather than and separately, and has the form Kern03 (); note1 ()


where . The first term in Eq. (1) is the hard-sphere (HS) contribution


while the second term can be factored into an isotropic square-well (SW) tail


modulated by an angle-dependent factor


The unit vectors are defined by the spherical coordinates in an arbitrarily oriented coordinate frame. Here we will put , where is Boltzmann’s constant and is the absolute temperature, and introduce the particle density . We use reduced units for temperature, , and density, , in the description of the thermodynamics. The above potential then ensures a proper bonding of the two particles depending upon the relative orientation and distance of the attractive caps on each sphere.

The square of the total coverage can be computed in terms of as


where is the Heaviside step function, equal to if and if , and where we have introduced the angular average


The integral can be readily evaluated to give Kern03 ()


Knowledge of the exact result (7) of integral (5) is then exploited to optimize the discretization of the angular integration appearing in all successive integral equations, illustrated in the next Section.

Iii Molecular integral equation approach

In the case of spherically symmetric potentials, the way to extract the thermophysical properties of a fluid has a long and venerable tradition in integral equation theory. Its central aim is the calculation of the pair distribution function , also typically computed in numerical simulations, from the pair potential . It is useful as well to introduce the total correlation function and the so-called direct correlation function defined through the Ornstein-Zernike (OZ) equation


An exact, albeit formal, relation holds between such functions and the pair potential:


where the last term in the argument of the exponential is a (non-explicit) functional of the correlation function, generally called a bridge function for historical reasons. Green () All the existing approximations may be recast into the form of an approximate bridge function in Eq. (9), the so-called closure equation. Most current algorithms also invoke the use of the auxiliary function , which is a continuous function even for discontinuous potentials such as hard spheres

The case of angle-dependent anisotropic potentials, although far more complex from an algorithmic point of view, follows essentially the same scheme. It was devised in the frame of molecular fluids Gray84 () and more recently adapted to the specific case of the Kern-Frenkel potential. Giacometti09b (); Giacometti10 () For completeness, the iterative procedure followed in Refs. Giacometti09b, ; Giacometti10, is briefly reviewed below.

iii.1 Iterative procedure

Our notation in this section will closely follow that of Gray and Gubbins in Ref. Gray84, , with only a prefactor difference; for instance, . Starting with a reasonable guess for the set of coefficients in the axial -frame, where , we use an expansion in spherical harmonics to obtain that in this frame depends only upon (),


where and the are spherical harmonics. Then we can use the closure relation


to obtain that, in this frame, still depends only upon (). The bridge function in this expression must be approximated, giving rise to such distinct closures as Percus-Yevick (PY) and hypernetted-chain (HNC); see below for the reference HNC (RHNC) closure used in this work. The inverse of an expansion like Eq. (10) is then used to compute the coefficients within the same frame,


To carry out Fourier transforms and so deconvolute the molecular OZ equation, Gray84 () we need to move at this point into an arbitrary space frame (often referred to as laboratory-frame) by means of a Clebsch-Gordan (CG) transform,


where the are Clebsch-Gordan coefficients. Fourier transforms then become Hankel transforms of the form


where is a spherical Bessel function of order . We can then return to a specific frame, the axial -frame, where this time . This can be achieved by means of an inverse Clebsch-Gordan transform,


Now one may use the Ornstein-Zernike equation in space, that in the axial -frame becomes


to obtain the new transform coefficients by matrix operations. As before, one needs now to return to a more general space frame through a Clebsch-Gordan transform in Fourier space,


because this allows the return to direct space by means of an inverse Hankel transform,


A final inverse Clebsch-Gordan transform then completes the return to the axial -frame we started with,


and thus yields a new estimate of the starting coefficients , in general different from the previous one. These steps are iterated until consistency between input and output coefficients is achieved. Table 1 summarizes the procedure.

Table 1: Schematic flow-chart for the solution of the OZ equation for the Kern-Frenkel angle-dependent potential. See Section III.1 for a description of the scheme.

iii.2 The RHNC closure and free energy

Although the second equation in this scheme, Eq. (11), is formally exact, it involves the calculation of the bridge function that in practice cannot be computed exactly, Hansen86 () as remarked earlier, and so an approximate closure is needed. Our approach is based on the RHNC approximation introduced in Ref. Lado73, for spherical potentials and later extended to molecular fluids. Lado82a (); Lado82b () Within this scheme, the closure equation takes on the assumed-known bridge function of a particular reference system to replace the actual unknown bridge function appearing in the exact closure. The goodness of the approximation clearly depends upon the quality of the chosen bridge function for the reference system. In the present case, for want of a better option, this is taken to be the hard-sphere model so that , where is the reference hard-sphere diameter. It has been demonstrated Rosenfeld79 (); Lado82 () that internal thermodynamic consistency can be improved upon treating as a variational parameter to be optimized. While the use of the hard-sphere bridge function is a natural assumption leading to a rather accurate approximation for spherically symmetric potentials, this is not as likely to be the case for a severely anisotropic potential such as the one-patch Kern-Frenkel model studied here. As we shall see below, this drawback is indeed confirmed by our findings, but better approximations for anisotropic potentials are not yet available.

Within the RHNC approximation, the excess free energy can be computed as Lado82b ()




In Eq. (22), is a Hermitian matrix with elements , , and is the unit matrix. In Eq. (23), directly expresses the RHNC approximation. Here is the reference system contribution, computed from the known free energy of the reference system as , with and calculated as above but with reference system quantities.

For the bridge function appearing in (23), we use the Verlet-Weis-Henderson-Grundke parametrization, Verlet72 (); Henderson75 () with the optimum hard sphere diameter selected according to a variational free energy minimization that yields the condition Lado82 ()


iii.3 Thermodynamics

The main strength of the RHNC closure hinges on the fact that, unlike most other closures, no further approximations are needed to obtain the free energy (as seen above) and other thermodynamic quantities. The pressure can be derived from a standard expression Gray84 () as


Introducing the cavity function and using the result


Eq. (25) becomes


which can be computed using Gaussian quadratures. Note that the second equality in Eq. (25) implies that the pressure depends upon the quality of , the other components being irrelevant.

The chemical potential can then be obtained from the exact thermodynamic relation


with the ideal quantities given by , , , where is the de Broglie wavelength.

Iv Improved Newton-Raphson algorithm

The iteration cycle described in Section III.1, wherein the output coefficients of one iteration directly become the input coefficients of the next, is known as Picard iteration. While obviously straightforward, it produces successive outputs that often converge only slowly or sometimes not at all, even for thermodynamic states that are known to exist. A standard remedy is to construct the new input coefficients for the next iteration as a damping linear combination of the current input and output sets. Broyles60 () We have implemented it in the efficient form proposed by Ng Ng74 () for generating a new input set of as an optimized linear superposition of the output sets from up to the previous four iterations.

But a more powerful procedure than such enhanced Picard cycles is available in the iterative application of Newton’s well-known root-finding algorithm. In the present context, however, Newton’s method, also known as the Newton-Raphson (NR) method, has the serious drawback of becoming so computationally intensive as to be prohibitive in practice, even for spherically symmetric models with just one coefficient. A clever meld of these two iteration techniques, producing a Newton-Raphson/Picard hybrid, was first proposed by Gillan Gillan79 () for spherically symmetric models, using a small number of so-called roof functions to represent the “coarse” features of for NR processing. (Here is the grid interval in the discrete space used in a numerical solution; the total number of grid points is .) Later, Labík, Malijevský, and Voňka (LMV) Labik85 (); Lomba89 () suggested an elegant alternative based instead on the NR processing of a small number, up to some cutoff , of values, where is the grid interval in space. In this work, we have implemented the LMV hybrid, but for just the coefficient, which makes the biggest contribution to , as explicitly illustrated by the results presented in Sec. VI.4, while the other coefficients are treated by a standard Picard cycle. Not only does the algebra become unwieldy if more components are included in the NR iterations, but for the Kern-Frenkel potential there is no obvious basis for choosing which additional components to include. We wish then to solve the one-component OZ equation (see Eq. (16))


for on the discrete grid, from to , where . Let be the desired solution, so that


and is a function of all the . If is our current value for the unknown, then we need to find the correction such that This is accomplished in the NR root-finding method by setting


Matrix inversion of Eq. (31) for the first points then produces the desired corrections .

Tracking the simplified one-component version of the Picard cycle in Section III.1, , leads to Labik85 ()


and completes the NR prescription. The discrete version lado67 () of the reciprocal Fourier transforms requires that the intervals and satisfy . In the present calculations, we have used , , and .

V Barker-Henderson thermodynamic perturbation theory

Barker-Henderson perturbation theory Barker67 (); Henderson71 (); Barker76 () hinges on the splitting of the Kern-Frenkel potential, Eq. (1), into the hard-sphere contribution, Eq. (2), and the remaining “perturbation” term,


This allows the high-temperature expansion of the free energy as


where is the free energy of the hard-sphere reference system, and where the first-order term,


can be easily computed in terms of the radial distribution function of the HS reference system; here is the hard-sphere packing fraction. The second-order term is, on the contrary, a highly non-trivial calculation involving higher-order correlation functions. An extension of the original Barker-Henderson alternative scheme yields the corresponding compressibility approximation that reads Gogelein12 ()


where is the reduced pressure of the HS reference system in the Carnahan-Starling approximation. CS69 () From here, pressure and chemical potential can be computed from the exact thermodynamic relations


Vi Results

vi.1 Pair distribution function

Unless otherwise stated, our results refer to , as in Ref. Giacometti09a, . Consider as initial state a reduced temperature for which the fluid is in a single phase at high density for all coverages examined here. We seek to determine the effect on the pair distribution function of reducing the coverage for the given state point. This is reported in Fig. 2 for three representative orientations: head-to-tail (HT), perpendicular (), and head-to-head (HH), corresponding to angles between the corresponding patch orientation vectors, respectively. (Similar plots were also considered in related systems, such as spherocylinders; see for instance Ref. Martinez-Haya03, )

Figure 2: The distribution function as a function of for three orientations of the patches: HT, ; , ; HH, , and different coverages from (square-well) to (Janus).

Clearly, while for the HT case is only mildly affected within the well, , both the and the HH pair distribution functions display a significant increase close to the contact point .

On the other hand, the coexistence lines progressively shift to lower temperatures for decreasing coverages, as we will see, and hence a fixed state point in the temperature-density plane is correspondingly moving relatively farther and farther from them, as coverage decreases.

In order to account for this and make different coverages comparable, we consider different state points that are comparably close to the gas-liquid coexistence lines. These are shown in Fig. 3 for decreasing coverage from to and two specific state points, side by side, that have different temperatures for the different coverages. In each case, we have first considered the largest computed density ( for all coverages) and the corresponding lowest computed temperature (decreasing with decreasing coverage). The panels on the left side of Fig. 3 correspond to state points expected to lie in the liquid phase at the respective coverages and are shown for decreasing coverage from top to bottom. The other set of chosen state points in the right-hand panels of Fig. 3 are all points lying in the respective gas phases (low temperatures and low densities) and are depicted again for decreasing coverage from top to bottom. In all cases, three different curves are reported corresponding to the HT, , and HH orientations of the two patches.

Figure 3: The distribution function as a function of for three orientations of the patches: HT, ; , ; HH, . In all cases, we present results for the highest and lowest densities studied at the lowest temperatures achieved at each coverage. From top to bottom, this corresponds to: , , (left), (right); , , (left), (right); , , (left), (right).

Consider first the high-density state points on the left. Few general features are readily apparent. In all cases, the HT curve exhibits a hard-spheres behavior with no discontinuity at the well edge, , as expected from the definition of the Kern-Frenkel potential.

Note that the value of this at contact, , decreases as the coverage decreases, since it becomes less and less likely to find particles with the HT orientation of the patches as decreases (further note the change in scale among different cases). Conversely, both and HH curves exhibit the usual discontinuity at , indicating that they are involved in bonding, with a progressive increase of the at contact, , as coverage decreases that is more marked in the HH than in the case.

A rather interesting pattern emerges from the low-density plots of the right-hand panels. Those are the cases where one expects an increase in micellization as coverage decreases. This is indeed confirmed by the results. As coverage decreases, the general trend is a significant increase of at contact, , the largest increase pertaining to the HH orientations, as expected. This clearly indicates the formation of clusters (micelles or vesicles) with an increasing fraction of saturated bonds. In particular, in the Janus case () the HT orientation gives a flat curve around , indicating an almost ideal behavior that reflects the almost complete absence of such orientations. However, we have observed no significant discontinuity on passing from to coverages that would indicate anomalous behavior of the Janus case. Therefore, RHNC is clearly not able to capture this effect with the present spherically-symmetric approximation of .

vi.2 Angular distributions

Complementary to previous cases, here we focus on the dependence of on just the orientations of and relative to within the square-well region. The expansion in spherical harmonics of in an arbitrary space frame reads


where we have introduced the rotational invariants Gray84 (); Hansen86 ()


In Ref. Giacometti10, , it was shown that upon defining


the resulting function of the polar coordinate of and the polar coordinate of the second patch reads


given that the axis is aligned with the patch of particle 1.

The behavior of as a function of is reported in Fig. 4 for three different orientations of the patches: HT (), (), HH (), and different coverages from to . The same high and low densities state points used before have been considered here. This identifies the preferential angular positions of the various different patch orientations.

Figure 4: The angular distribution as a function of for three orientations of patch (, , ) given that patch is pointing up (): In all cases, we present results for the highest and lowest densities studied at the lowest temperatures achieved at each coverage. From top to bottom, this corresponds to: , , (left), (right); , , (left), (right); , , (left), (right). The circular arrowed insets refer to the patch orientation, with (always up) and rotating.

Consider the high density state point first, depicted in the left-hand panels of Fig. 4 for decreasing coverages from top to bottom. State points are the same discussed in Fig. 3. For sufficiently large patches (, not shown here), the only significant peak in the distribution is observed for and . For such high coverages, HH alignments are uniformly distributed along all solid angles (remember that there is azimuthal symmetry), whereas HT alignment is preferentially found in the backward direction, .

The situation changes as the coverage decreases from , with the development of further peaks for perpendicular orientation of the patches () at and for head-to-head orientation of the patches () at . The physical interpretation of these results is that, under high density and low temperature conditions, head-to-tail (HT) and head-to-head (HH) alignments of the patches are preferentially found for particles in the transversal direction, , for low coverages ().

Next, we consider the low density points reported in the right-hand panels of Fig. 4, again for decreasing coverages from top to bottom. Unlike the previous case, we find a clear predominance of the HH antiparallel alignment in the forward direction () and modulated layering for both HT and patch orientations that become increasingly structured as coverage decreases. These results can be contrasted with the analogous results given in Ref. Giacometti10, for the two-patch case and extend those given there for only high and low coverages. The layering is a clear reflection of an increasing tendency to micellization, in agreement with numerical simulation results.

vi.3 Coefficients of rotational invariants

In this section, we follow the the notations already introduced in our previous work. Giacometti10 () The coefficients of rotational invariants are


where the are rotational invariants. Here we have explicitly considered the first coefficients occurring in the multipole expansion Gray84 () that account up to quadrupole-quadrupole interactions. Stell81 ()

Explicit expressions for the first few are Stell81 ()


Other expressions can be found in Ref. Stell81, . We note that used in past work Giacometti10 () and further that .

In Appendix A, we explicitly derive Eq. (45) for two specific and representative cases. Some of the coefficients have particularly interesting physical interpretations: the term is the coefficient of ferroelectric correlation, the term the coefficient of dipolar correlation, the term the coefficient of nematic correlation, and so on.

The results for these coefficients are reported in Fig. 5, with the same ordering as before. Hence the left-hand panels show plots of the high-density points and decreasing coverage, while the right-hand panels depict plots of the low-density points and again decreasing coverages. Plots on the same side have been drawn to the same scale so that differences may be readily appreciated.

The high-density plots (left-hand panels) have hardly any dependence on the particular projection, as could have been guessed from the outset. With , we clearly find correlations (that is, non-vanishing coefficients) only within the well, , along with and negatively correlated, positively correlated, and almost uncorrelated. Similar behavior occurs for the low-density state points where, however, the correlation within the well is approximately constant, with , and . Note that in the last, Janus case (), the and ordering appear to be inverted, signaling an incomplete agreement with the other cases, likely due to an insufficient lowering of the temperature, in agreement with previous findings of Sections VI.1.

Figure 5: The rotational invariants as a function of for several triplets. In all cases, we present results for the highest and lowest densities studied at the lowest temperatures achieved at each coverage. From top to bottom, this corresponds to: , , (left), (right); , , (left), (right); , , (left), (right), as in Fig. 3.

Next we consider a second set of coefficients given by , , , , . These are reported in Fig. 6 with the same distribution as before. Even in this case, all coefficients have non-vanishing values within the well and have thus been plotted to the same scale. Again, the trend appears to be rather clear, with the coefficient negative with decreasing contact values for decreasing patch size, indicating an increasing anticorrelation in the respective orientations as coverage decreases; also has negative value, whereas all others coefficients present positive values indicating positive correlations. This is true for both high- and low-density states.

Figure 6: The rotational invariants as a function of for several triplets. In all cases, we present results for the highest and lowest densities studied at the lowest temperatures achieved at each coverage. From top to bottom, this corresponds to: , , (left), (right); , , (left), (right); , , (left), (right).

vi.4 RHNC molecular reference angular components of radial distribution functions and MC results

In order to assess the structural results previously discussed, in this section we compare directly the RHNC molecular-frame spherical harmonic coefficients and MC results.

Figure 7: The