Free-Energy Functional Method for Inverse Problem of Self Assembly
aa r X i v : . [ c ond - m a t . s o f t ] A p r Free-Energy Functional Method for Inverse Problem of Self Assembly
Masashi Torikai a) Department of Physics Engineering, Faculty of Engineering, Mie University, 1577 Kurimamachiya-cho, Tsu,Mie 514-8507, Japan (Dated: 9 April 2015)
A new theoretical approach is described for the inverse self-assembly problem, i.e., the reconstruction of theinterparticle interaction from a given structure. This theory is based on the variational principle for thefunctional that is constructed from a free energy functional in combination with Percus’s approach [J. Percus,Phys. Rev. Lett. , 462 (1962)]. In this theory, the interparticle interaction potential for the given structureis obtained as the function that maximizes the functional. As test cases, the interparticle potentials fortwo-dimensional crystals, such as square, honeycomb, and kagome lattices, are predicted by this theory. Theformation of each target lattice from an initial random particle configuration in Monte Carlo simulations withthe predicted interparticle interaction indicates that the theory is successfully applied to the test cases. I. INTRODUCTION
In condensed phases, diverse interactions between theconstituents (such as atoms, molecules, micelles, nano-and microparticles) produce a variety of simple andcomplex structures. Besides a plethora of the crystal,liquid crystal, and quasicrystal structures formed by con-ventional atoms and molecules, square and honeycomblattices of colloidal nanocrystals, kagome lattice of tri-block Janus particles, and quasicrystals of dendrimers and binary nanocrystals are additional experimentalexamples of non-trivial structures. Molecular simulationsfor systems with several interparticle interactions havebeen performed to enumerate the stable structures; it hasbeen found that particles with a rather simple interactionpotential can assemble into complex structures, e.g.,hard sphere (HS) plus linear ramp, HS plus squarewell, HS plus square shoulder, and Lennard-Jonesplus Gaussian (LJG) potential particles assemble intoa two-dimensional (2D) quasicrystal. Density functionaltheory (DFT), another tool to find the thermodynami-cally stable phases, has shown that in three dimensions(3D), the LJG fluid forms a body centered cubic (bcc)crystal as a stable phase. When a structure of a condensed phase is known fromexperiments or is designed artificially, determination ofthe interparticle interaction is more efficient by using theinverse approach, in which the interparticle interactionis derived from the given structure, than using theforward approach of exhaustive search by the molecularsimulations or DFT. The inverse statistical mechanicalmethod, which is the most elaborated inverse approach,has been successfully used to design the interparticleinteraction potentials that generate the square, hon-eycomb, kagome, and rectangular lattices in2D, and the simple cubic, bcc, simple hexagonal, diamond, and wurtzite lattices in 3D. In theinverse statistical mechanical method, the interparticleinteraction is optimized with respect to the energy and a) [email protected] mechanical stability criteria via molecular simulations.The studies on the inverse problem at finite temperaturehave mostly been based on the molecular simulationsat finite temperature. (It was suggested in Ref. 18that the molecular simulations at finite temperature maynot be necessary. This is based on the observationthat the target structure shows good stability at finitetemperature if the interaction potential is optimizedat zero temperature via the minimization of a specificsimulated annealing energy defined in Ref. 16. However,this strategy is also dependent on computer simulations.)For the liquid target phases, the reverse Monte Carlo(MC) method, which also uses the MC simulation inthe optimization, is another successful inverse approach.Here the aim is to formulate a theory to reconstruct theinterparticle potential that can stabilize a given targetstructure (more specifically, the single- and two-particledistributions) at finite temperature. In this paper, Ipresent a new simulation-free inverse method that is avariational method based on the free-energy functionaltheory akin to DFT. The interparticle interaction func-tion is defined as the function that gives the maximumof the functional. I applied this method to the square,honeycomb, and kagome lattices, and obtained the corre-sponding interparticle potentials. The potentials are thenused in a series of simulated annealing MC calculationsstarting from random configurations. In most cases, theresulting solid contains a few grain boundaries and manydefects. However, it is found that for each predictedpotential, in at least a few percent of the simulationsthe particles spontaneously form the appropriate targetlattice with a small number of defects. Although thesmall success rate in the simulations indicates that theinteraction potentials obtained here are not entirely opti-mized as the potentials obtained by the inverse statisticalmechanical method, the observed self-assembly into thetarget lattice implies that the method introduced hereis another promising approach to the inverse problem ofthe self-assembly. II. THEORY
We consider a single component system in a d -dimensional space. The system comprises particles in-teracting with a pairwise-additive potential v ( x ). Thegrand potential of the system, with the temperature T ,chemical potential µ , and volume V , in the presence ofan external field φ ext ( x ), is defined asΩ[ ϕ ] = − β − ln Ξ[ ϕ ] , (1)where β = 1 /k B T is the inverse temperature, ϕ ( x ) = µ − φ ext ( x ) is the intrinsic chemical potential, and Ξ[ ϕ ]is the grand partition functionΞ[ ϕ ] = ∞ X N =0 Z d r ( N ) λ dN N ! × exp (cid:20) − β N − X i N X j>i v ( r i − r j ) + β N X i ϕ ( r i ) (cid:21) , (2)where r i and λ denote the coordinate of the i th particleand the de Broglie thermal wavelength, respectively.The key functional for this approach is˜ A [ ρ, ψ ] = Ω[ ψ ] + Z ρ ( x , [ ϕ ]) ψ ( x )d x , (3)where ρ ( x , [ ϕ ]) denotes the single-particle density in thepresence of the intrinsic chemical potential ϕ ( x ) and ψ ( x ) is an independent function. The maximum of˜ A [ ρ, ψ ] with respect to ψ , i.e., A [ ρ ] = sup ψ ˜ A [ ρ, ψ ] , (4)is the intrinsic free energy, the Legendre transformationof the grand potential Ω[ ψ ]. The maximum of ˜ A [ ρ, ψ ] isachieved when ψ is equivalent to the intrinsic chemicalpotential ϕ , i.e., A [ ρ ] = ˜ A [ ρ, ϕ ]. The inequality relation˜ A [ ρ, ϕ ] ≥ ˜ A [ ρ, ψ ] for any ψ ( x ) (5)provides us the variational method to determine theintrinsic chemical potential ϕ ( x ) that gives rise to agiven density profile ρ ( x , [ ϕ ]): the ψ ( x ) that maximizesthe functional ˜ A [ ρ, ψ ] is the intrinsic chemical potential ϕ ( x ). The variational method for ˜ A [ ρ, ψ ], in which theproper ϕ ( x ) is obtained for a given ρ ( x ), can be viewedas the inverse of the DFT, in which the proper ρ ( x )is obtained for a given ϕ ( x ).The variational method for ˜ A [ ρ, ψ ] can be adapted toobtain the information about the interaction potential v ( x ) using Percus’s idea. If a particle is fixed at theorigin of the coordinate system, the remaining particlesfeel the external field v ( x ) due to the fixed particle inaddition to the original external field φ ext ( x ). In thissituation, the single-particle density becomes ρ ( x , [ ϕ fix ]),where ϕ fix ( x ) = µ − [ φ ext ( x ) + v ( x )] = ϕ ( x ) − v ( x ) is the intrinsic chemical potential in the presence of the fixedparticle. Percus showed that ρ ( x , [ ϕ fix ]) is related tothe single- and two-particle density in the absence of thefixed particle by ρ ( x , [ ϕ fix ]) = ρ (2) ( x , , [ ϕ ]) ρ (0 , [ ϕ ]) , (6)where ρ (2) ( x , x ′ , [ ϕ ]) is the two-particle density between x and x ′ in the absence of the fixed particle. As discussedabove, the function ψ that maximizes ˜ A [ ρ ( x , [ ϕ fix ]) , ψ ] is ϕ fix ( x ). Combining the variational method for ˜ A [ ρ, ψ ]with Percus’s idea, we find that the function ψ ( x ) thatmaximizes ˜ A (cid:2) ρ (2) (cid:14) ρ, ψ (cid:3) is ϕ fix ( x ). Thus, for a given setof ρ ( x ) and ρ (2) ( x , x ′ ), we can obtain the interaction po-tential v ( x ) = ϕ ( x ) − ϕ fix ( x ) through the maximizationof ˜ A (cid:2) ρ (2) (cid:14) ρ, ψ (cid:3) with respect to ψ ( x ).The inhomogeneous version of Percus’s relation (6)can be derived along the same line as the derivation ofthe homogeneous version in Ref. 22 by simply removingthe assumption of the homogeneity. In an inhomoge-neous and symmetry-broken phase, however, the naivedefinition of the distribution functions loses its unique-ness. In some systems, for example, a ferromagnet,the uniqueness of the order parameter is recovered byapplying a field that breaks the symmetry of the originalsystem, and the order parameter without the field can beobtained by taking zero-limit of the symmetry-breakingfield followed by the thermodynamic limit. I do not knowwhether the same procedure applies to the distributionfunctions of inhomogeneous fluids, thus cannot providea rigorous proof of Percus’s relation for inhomogeneouscases; however, I postulate the validity of Eq. (6) in thispaper. III. APPLICATIONS TO 2D CRYSTALS
The free-energy functional method introduced in Sec.II is applicable to any 2D and 3D system that ischaracterized by its single- and two-particle densities.In this paper, as test cases, I use 2D crystals suchas square, honeycomb, and kagome lattices as targetstructures. While the interaction potentials that stabilizethese target crystals have already been found in previousstudies and are currently of little novelty; never-theless, these crystals still serve as good test cases.From the variational principle for ˜ A (cid:2) ρ (2) (cid:14) ρ, ψ (cid:3) , itfollows that the intrinsic chemical potential ϕ fix ( x )is the solution of the Euler-Lagrange equation δ ˜ A (cid:2) ρ (2) (cid:14) ρ, ψ (cid:3) /δψ ( x ) = 0, i.e., ρ (2) ( x , x ′ ) ρ ( x ′ ) = − δ Ω (cid:2) ψ (cid:3) δψ ( x ) (cid:12)(cid:12)(cid:12)(cid:12) ψ = ϕ fix . (7)A direct solution of this equation is computationallydemanding. Therefore, in this paper, I use a trialfunction v tr ( x , { x i } ) with a set of variational parameters { x i } as the interaction potential. The function ψ ( x )then becomes ψ tr ( x , { x i } ) = µ − v tr ( x , { x i } ), and thefunctional ˜ A [ ρ (2) (cid:14) ρ, ψ (cid:3) is reduced to the function of { x i } .The set of parameters { x i } that maximizes ˜ A (cid:2) ρ (2) (cid:14) ρ, ψ tr (cid:3) defines the required interaction potential.In the absence of the external field, the equilibriumconfiguration is determined solely by βv ( x ), becausethe grand partition function (2) depends on v ( x ) onlythrough βv ( x ). Thus, the optimal interparticle potentialis linearly dependent on the temperature.The grand potential Ω[ ψ tr ] can be expanded in afunctional Taylor expansion. Since ψ tr ( x ) for the typ-ical interparticle interactions diverges around x = 0due to the short range repulsion, the activity z ( x ) =exp[ βψ tr ( x )] /λ is easier to handle than ψ tr ( x ) itself.The functional Taylor expansion of the grand potentialin powers of the deviation of activity ∆ z ( x ) = z ( x ) − z is β Ω[ ψ tr ] = − ∞ X n =0 n ! Z · · · Z δ n ln Ξ[ ϕ ] δz ( x ) δz ( x ) . . . δz ( x n ) (cid:12)(cid:12)(cid:12)(cid:12) z = z n Y i =1 ∆ z ( x i )d x i , (8)where z = e βµ /λ is the activity in the absence ofthe fixed particle. The functional derivatives of ln Ξwith respect to z ( x ) can be expressed in terms of themultiparticle distribution functions. Substituting (6)and (8) into (3), we obtain β ˜ A (cid:2) ρ (2) (cid:14) ρ, ψ tr (cid:3) = β Ω[ ϕ ] + β Z ρ (2) ( x , ρ (0) ψ tr ( x )d x − Z ρ ( x ) ζ ( x )d x − Z Z (cid:2) ρ (2) ( x , x ′ ) − ρ ( x ) ρ ( x ′ ) (cid:3) ζ ( x ) ζ ( x ′ )d x d x ′ , (9)up to the second order in ∆ z ( x ), where ζ ( x ) =∆ z ( x ) /z = exp[ − βv tr ( x )] −
1. The approximationof the truncated Taylor series is justified if ζ ( x ) issmall. Therefore, when v tr ( x ) has a negative value, thisapproximation may break down at very low temperature.Hereafter, I restrict the consideration to a case in whichthe minimum of βv tr ( x ) is − v tr ( x ), I use the LJG potential,which is an isotropic potential constructed by adding aGaussian function to the Lennard-Jones potential. Here,the trial function is the LJG potential with positive andnegative Gaussians: v tr ( x ) = ǫ ((cid:18) αx (cid:19) − (cid:18) αx (cid:19) − c exp (cid:20) − ( αx − x ) σ v (cid:21) + c exp (cid:20) − ( αx − x ) σ v (cid:21)(cid:27) , (10)where σ v , c and c denote the Gaussian functions’ width,depth, and height, respectively. In this paper, I use the parameters σ v = √ . c = 1, and c = 4. Thevariational parameters x i ( i = 1 ,
2) are used to adjustthe positions of the Gaussian functions. The valuesof the remaining parameters α and ǫ are determinednumerically so that the first minimum of v tr ( x ) will be atthe nearest neighbor atomic distance and the value of theglobal minimum will be −
1. Considering these criteria forthe trial function, the functional ˜ A [ ρ (2) (cid:14) ρ, ψ tr (cid:3) is reducedto a function ˜ A ( x , x ).The atomic positions in the target perfect crystals, { a i } , are the input for the free-energy functional method.I use units for which the nearest neighbor atomic distanceis unity. In the present paper, I assume that the i thparticle fluctuates around a i and is independent of theother particles, and that its fluctuation is given bya Gaussian distribution. The single- and two-particledistribution functions are then ρ ( x ) = X i f ( x − a i ) , (11) ρ (2) ( x , x ′ ) = X i f ( x − a i ) X j = i f ( x ′ − a j ) , (12)where f ( x ) denotes the Gaussian distribution f ( x ) = 1 πσ exp (cid:20) − (cid:16) x σ (cid:17) (cid:21) . (13)The resulting interaction potential depends on the stan-dard deviation σ . In principle, any value of σ ispermissible in this method. However, because ρ ( x ) and ρ (2) ( x , x ′ ) with sharp peaks are not preferable for theaccuracy of the numerical integration in (9), the broadestbut physically acceptable distribution f ( x ) is better.Therefore, in this paper, σ is set to 0 .
15, which is thetypical value of the Lindemann ratio at crystallization.Since ρ (2) ( x , /ρ (0) remains finite near x = 0, where ϕ fix ( x ) diverges rapidly, the approximations (11) and(12) cause the second term at the right-hand side of (9)to diverge. This is caused by neglecting the correlationbetween the particles in the derivation of (12). Toinclude the non-negligible particle correlation due tothe strong repulsion between the particles separated byshort distances, I modified (12) by multiplying it byexp (cid:2) − β (cid:0) / | x | − (cid:1)(cid:3) for | x | <
1. The modification factoris proportional to the ideal gas density in a repulsiveexternal field 1 / | x | .One of the particles in the perfect crystal, for example a , is set at the origin of the coordinate system. Ingeneral, ρ (2) ( x ,
0) depends on the choice of a since eachparticle is not necessarily equivalent in the configurationof the target structure. In such a case, the averagevalues of ˜ A over all possible choices of a should beused. However, this is not necessary for the target crystalstructures used in this study.I have numerically determined ˜ A ( x , x ) according to(9) using MC integration as implemented in the VEGASalgorithm in GNU Scientific Library. (cid:82)(cid:88)(cid:121) (cid:82)(cid:88)(cid:107) (cid:82)(cid:88)(cid:57) (cid:82)(cid:88)(cid:101) (cid:82)(cid:88)(cid:51) (cid:107)(cid:88)(cid:121) (cid:107)(cid:88)(cid:107)(cid:82)(cid:88)(cid:121)(cid:82)(cid:88)(cid:107)(cid:82)(cid:88)(cid:57)(cid:82)(cid:88)(cid:101)(cid:82)(cid:88)(cid:51)(cid:107)(cid:88)(cid:121) x x − − − −
10 0
FIG. 1. The density plot of ˜ A ( x , x ) for the square lattice.The location of the peak in ˜ A ( x , x ) is first estimated on acoarse grid (∆ x i = 0 .
1) and is then estimated on a finer grid(∆ x i = 0 .
02) around the first estimated peak.TABLE I. The radial distances and coordination numbers forthe first four nearest neighbors for perfect lattices.Square Honeycomb Kagome TriangularFirst 1, 4 1, 3 1, 4 1, 6Second √
2, 4 √
3, 6 √
3, 4 √
3, 6Third 2, 4 2, 3 2, 6 2, 6Fourth √
5, 8 √
7, 6 √
7, 8 √
7, 12
The density plots of ˜ A ( x , x ) for target lattices areshown in Figs. 1, 2, and 3. The peaks in ˜ A ( x , x )for square, honeycomb, and kagome lattices are locatedat ( x , x ) = (2 . , . . , . . , . x ≈ . (cid:82)(cid:88)(cid:121) (cid:82)(cid:88)(cid:107) (cid:82)(cid:88)(cid:57) (cid:82)(cid:88)(cid:101) (cid:82)(cid:88)(cid:51) (cid:107)(cid:88)(cid:121) (cid:107)(cid:88)(cid:107)(cid:82)(cid:88)(cid:121)(cid:82)(cid:88)(cid:107)(cid:82)(cid:88)(cid:57)(cid:82)(cid:88)(cid:101)(cid:82)(cid:88)(cid:51)(cid:107)(cid:88)(cid:121) x x − − FIG. 2. The density plot of ˜ A ( x , x ) for the honeycomblattice. The estimation grid is the same as in Fig. 1. (third) nearest neighbor distance, because the secondminimum of the honeycomb (kagome) potential has alarger coordination number at the second (third) nearestneighbor distance than the third (second) one.The square lattice potential found in Ref. 11 has twominima at the first and fourth nearest neighbor distance,like the square lattice potential found here. But thesepotentials are not similar, because the former has a deepminimum at the first nearest neighbor distance and avery shallow minimum at the fourth nearest neighbordistance. The potential for honeycomb lattice here andthat in Ref.10 have some common features: the firstminimum has positive value and the second minimumis at the second nearest neighbor distance. The potentialfor the kagome lattice here is completely different fromthe known potentials for the kagome lattice, which arepurely repulsive. To examine the validity of these predicted potentials,I performed MC simulations for a constant number ofparticles N , volume (area) V , and temperature T . Foreach simulation, N was greater than 500. The simulationbox was a square of side L with periodic boundaryconditions; L = √ V was defined so that the density ofthe system N/V corresponds to that of the target perfectlattice. If N is not appropriate to fill the simulationbox with the unit cells of the target lattice, the successrate of the crystallization into the target lattice willdecrease. For square lattice potential, N was chosen tobe a “magic number,” with which the unit cells fill thesimulation box when one of the sides of the square unitcells and that of the simulation box are parallel. For the (cid:82)(cid:88)(cid:121) (cid:82)(cid:88)(cid:107) (cid:82)(cid:88)(cid:57) (cid:82)(cid:88)(cid:101) (cid:82)(cid:88)(cid:51) (cid:107)(cid:88)(cid:121) (cid:107)(cid:88)(cid:107)(cid:82)(cid:88)(cid:121)(cid:82)(cid:88)(cid:107)(cid:82)(cid:88)(cid:57)(cid:82)(cid:88)(cid:101)(cid:82)(cid:88)(cid:51)(cid:107)(cid:88)(cid:121) x x − − − FIG. 3. The density plot of ˜ A ( x , x ) for the kagome lattice.The estimation grid is the same as in Fig. 1. − v ( x ) x FIG. 4. The predicted LJG potentials for square (dashedblack curve), honeycomb (solid red curve), and kagome(dotted green curve) lattices. These are given by Eq. (10)with ( x , x ) = (2 . , . . , . . , . honeycomb and kagome lattice potentials, N was chosenso that the distortion of the lattice remains small whenthe target lattice fits in to the simulation box under theperiodic boundary conditions. A random configurationwas used as an initial configuration for each series ofcooling process. Ten simulation runs were performed forsquare lattice potential, and all resulting solid phases aresquare lattice with some defects and no grain bound-aries. For the honeycomb and kagome lattice potentials, FIG. 5. Results of the 529-particle MC simulations for theLJG potential with ( x , x ) = (2 . , . k B T = 2 . k B T = 1 . N/V = 1 in the square simulationbox of linear dimension L = 23 .
0. The particle pairs separatedby a distance smaller than 1 . most simulation runs resulted in solid phases with grainboundaries and many defects; however, these solid phasesexhibited local structures that resembled those of thetarget lattice. Several runs (six and four runs, out ofthirty runs, for honeycomb and kagome lattice potential,respectively) resulted in crystalline structures withoutgrain boundaries and only a few defects. The snapshotsof these phases are shown in Figs. 5, 6, and 7, and thefigures clearly show that the particles interacting with theLJG potentials predicted by the free-energy functionaltheory self-assemble into the target crystals. IV. DISCUSSION AND SUMMARY
In this paper, I developed the free-energy functionalmethod for the reconstruction of the interparticle inter-action potential from a given structure by combiningthe variational principle for the functional ˜ A [ ρ, ψ ] withPercus’s idea. In this free-energy functional method, thedesired interaction potential is given as the function thatmaximizes ˜ A (cid:2) ρ (2) (cid:14) ρ, ψ (cid:3) . This method was successfullyapplied to the square, honeycomb, and kagome lattices.The free-energy functional method introduced hererequires single- and two-particle densities, ρ ( x ) and ρ (2) ( x , x ′ ), as input. To obtain the interaction potentialcorresponding to an artificially designed structure, e.g.,the sets of atomic positions { a i } , ρ ( x ) and ρ (2) ( x , x ′ )must be designed or approximated as was done in this FIG. 6. Results of the 544-particle MC simulations for theLJG potential with ( x , x ) = (1 . , . k B T = 2 . k B T = 0 .
40 at
N/V = 4 / √ L = 26 . work. If, instead, the method is used to obtain a modelpotential that reproduces an experimentally observedstructure, experimental data for ρ ( x ) and ρ (2) ( x , x ′ ) aresufficient.The self-assembly of the target lattices in the MCsimulations shows that the free-energy functional methodworks properly. However, this does not necessarilymean that the interparticle interactions obtained hereare entirely optimized for the assembly of the targetlattices. Indeed, if we perform several MC simulationsusing ( x , x ) = (1 . , . A in Eq. (9), the use of the Gaussian approximationfor ρ ( x ) in Eq. (11), and the independent-fluctuationapproximation for ρ (2) ( x , x ′ ) in Eq. (12). More so-phisticated approximations will improve the resultinginterparticle interaction potentials. The maximization of˜ A in ( x , x ) space may also give rise to the unsatisfactoryquality of the resulting structures. The maximization onthe full functional space of interparticle interactions isnecessary for the best prediction within the frameworkof the present theory.The Taylor expansion of ˜ A not only affects the accu-racy of the prediction but also the range of applicabilityof this method. As discussed in Sec. III, the choice of thetemperature and the trial potential is restricted to justify FIG. 7. Results of the 504-particle MC simulations for theLJG potential with ( x , x ) = (1 . , . k B T = 2 . k B T = 0 . N/V = 3 / √ L = 24 . the truncation of the Taylor series. Therefore, the presentform of this method does not provide a unified descriptionof the optimal interaction potential for a given structureover a wide range (including T = 0) of temperature.Unfortunately, better approximations for ˜ A [ ρ, ψ ] otherthan the Taylor series are not available yet.Even within the truncated Taylor series approxima-tion, the restriction on the temperature can be relaxedusing non-negative trial potentials. Although this is anadditional strong restriction on the interaction potential,we can expect that non-negative interaction potentialssuffice for a wide variety of target structures at finitetemperature, considering that such potentials producevarious ground state structures. The interaction potential considered here is restrictedto the isotropic one for simplicity. It was found tobe sufficient for the self-assembly of the target crystalsconsidered here, as expected from the past work for theground state.
The class of target crystals that cannotbe assembled by isotropic particles is not yet clear, butthe isotropy-assumption of the potential will break downand must be discarded if the target structure is stronglyanisotropic.Improving the numerical efficiency is necessary forthe future work, e.g., the use of the trial potentialwith many variational parameters; the inverse problemfor the 3D structure, which requires time-consumingsix-dimensional numerical integration in Eq. (9). Theuse of a terraced (discretized) interparticle potential is apossible candidate for the way to reduce the computationtime. The terraced potentials are used in molecularsimulations and analytical calculations as anefficient and satisfactorily accurate way to investigate themany body systems.The theory introduced here is also applicable to theinverse problem for liquids; in this case, the interparticleinteraction is reconstructed from a given pair distributionfunction g ( x , x ′ ) or a radial distribution function g ( r ).In fact, the method is more suitable for homogeneoussimple liquids than for crystals, because for liquids, theapproximations (11) and (12) for distribution functionsare unnecessary. This is because ρ ( x ) = ρ is constant and ρ (2) ( x , x ′ ) is equal to ρ g ( x , x ′ ) in homogeneous simpleliquids. While the inverse problem for liquids has beensolved by the reverse MC method, our free-energyfunctional method serves as a new theoretical approachto tackle this problem. ACKNOWLEDGMENTS
I would like to thank Professor Akira Yoshimori forhelpful discussions. W. H. Evers, B. Goris, S. Bals, M. Casavola, J. de Graaf, R. vanRoij, M. Dijkstra, and D. Vanmaekelbergh, Nano letters ,2317 (2013). Q. Chen, S. C. Bae, and S. Granick, Nature , 381 (2011). X. Zeng, G. Ungar, Y. Liu, V. Percec, A. E. Dulcey, and J. K.Hobbs, Nature , 157 (2004). D. V. Talapin, E. V. Shevchenko, M. I. Bodnarchuk, X. Ye,J. Chen, and C. B. Murray, Nature , 964 (2009). E. Jagla, Physical Review E , 1478 (1998). A. Skibinsky, S. Buldyrev, A. Scala, S. Havlin, and H. Stanley,Physical Review E , 2664 (1999). T. Dotera, T. Oshiro, and P. Ziherl, Nature , 208 (2014). M. Engel and H.-R. Trebin, Physical Review Letters , 225505(2007). A. Suematsu, A. Yoshimori, M. Saiki, J. Matsui, and T. Odagaki,The Journal of Chemical Physics , 244501 (2014). M. Rechtsman, F. Stillinger, and S. Torquato, Physical ReviewLetters , 228301 (2005). M. Rechtsman, F. Stillinger, and S. Torquato, Physical ReviewE , 011406 (2006). E. Marcotte, F. H. Stillinger, and S. Torquato, The Journal ofChemical Physics , 164105 (2011). A. Jain, J. R. Errington, and T. M. Truskett, Physical ReviewX , 031049 (2014). G. Zhang, F. H. Stillinger, and S. Torquato, Physical Review E , 042309 (2013). M. Rechtsman, F. Stillinger, and S. Torquato, Physical ReviewE , 021404 (2006). A. Jain, J. R. Errington, and T. M. Truskett, Soft Matter ,3866 (2013). M. Rechtsman, F. Stillinger, and S. Torquato, Physical ReviewE , 031403 (2007). A. Jain, J. R. Errington, and T. M. Truskett, The Journal ofChemical Physics , 141102 (2013). R. L. McGreevy and L. Pusztai, Molecular Simulation , 359(1988). A. Lyubartsev and A. Laaksonen, Physical Review E , 3730(1995). J.-M. Caillol, Journal of Physics A: Mathematical and General , 4189 (2002). J.-P. Hansen and I. R. McDonald,
Theory of Simple Liquids , 4thed. (Academic Press, Oxford, 2013). J. Percus, Physical Review Letters , 462 (1962). M. Galassi, J. Davies, J. Theiler, B. Gough, G. Jungman,P. Alken, M. Booth, F. Rossi, and R. Ulerich,
GNU ScientificLibrary Reference Manual , 3rd ed. (Network Theory Limited.,2013). H. Cohn and A. Kumar, Proceedings of the National Academyof Sciences of the United States of America , 9570 (2009). G. A. Chapela, L. E. Scriven, and H. T. Davis, J. Chem. Phys. , 4307 (1989). L. Xu, S. V. Buldyrev, C. A. Angell, and H. E. Stanley, PhysicalReview E , 031108 (2006). M. N. Bannerman, R. Sargant, and L. Lue, Journal of Compu-tational Chemistry , 3329 (2011). K. B. Hollingshead, A. Jain, and T. M. Truskett, Journal ofChemical Physics , 1 (2013).30