Calculating Binodals and Interfacial Tension of Phase-Separated Condensates from Molecular Simulations with Finite-Size Corrections

Abstract

We illustrate three methods for calculating the binodals of phase-separated condensates from molecular simulations. Because molecular simulations can only be carried out for small system sizes, correction for finite sizes may be required for making direct comparison between calculated results and experimental data. We first summarize the three methods and then present detailed implementation of each method on a Lennard-Jones fluid. In the first method, chemical potentials are calculated over a range of particle densities in canonical-ensemble simulations; the densities of the dilute and dense phases at the given temperature are then found by a Maxwell equal-area construction. In Gibbs-ensemble Monte Carlo, the exchange between separated dilute and dense phases is simulated to obtain their densities. Lastly, slab-geometry molecular dynamics simulations model the dilute and dense phases in coexistence and yield not only their densities but also their interfacial tension. The three types of simulations are carried out for a range of system sizes, and the results are scaled to generate the binodals corrected for finite system sizes. Size-corrected interfacial tension is also produced from slab-geometry molecular dynamics simulations.

Keywords: Binodals; Biomolecular condensates; FMAP; Finite-size scaling; Interfacial tension; Lennard-Jones fluid; Molecular dynamics simulation; Monte-Carlo simulation; Phase separation.

Fig. 1

The Lennard-Jones (LJ) potential and its shifted version, for the interaction between a pair of particles. The inset shows a zoom into a small range of the interparticle distance around the cutoff r_cut = 3, showing that LJ shifted is zero at r = r_cut. For r > r_cut, LJ shifted is fixed at 0. Note that, at r = r_cut, although LJ shifted is continuous, its slope and hence the interparticle force have a discontinuity.

Fig. 2.

Illustration of calculating the chemical potential by Widom insertion. (a) The inserted particle has interaction energies u_I1, u_I2, …, and u_IN with particles 1, 2, …, and N in the system under constant N, V, and T. Note that the insertion is fictitious because the inserted particle merely probes its interactions with the N particles but does not affect their motions. (b) Monte Carlo simulations of the N-particle system. A particle is randomly selected from the list of N, and a displacement is then proposed, with uniform probability within a cube centered at the original position and with side length 2l_max (dashed box). The would-be change in the total energy, ΔU, of the system is calculated, and the proposed displacement is accepted with the probability shown.

Fig. 3.

Illustration of FMAP for a two-dimensional system. (a) Top: the “potential” function f(x) of the inserted particle, shown in a color scale with red and blue representing positive and negative values, respectively. Bottom: the “charge” density g(x) of the N-particle system. The black dot displays the position of a particle, and the light red shading represents the smear of a point charge. (b) Top: the discretized version of panel (a) top, on a 16 × 16 grid, with a grid spacing of 0.5. Bottom: the discretized version of panel (b) bottom. The point charge of each particle is distributed to the four nearest pixels. (c) Top: after the Fourier transform (FT) of f(x) and g(x), the product F*(k)G(k) undergoes inverse Fourier transform (IFT) to yield the energy function U_I(x_I), for the interaction of the inserted particle at position x_I with the N-particle system. Bottom: with an increase in grid points from 16 × 16 to 80 × 80, the energy function becomes very accurate.

Fig. 4.

The excess chemical potential, calculated by FMAP, over a range of particle densities. (a) Fifteen systems, with increasing numbers of particles in the same volume, are each simulated at constant N, V, T. (b) The excess chemical potential at 15 densities, and its fit to a 5th-order polynomial. The numbers of particles are N_i = 30i, where i = 1 to 15; the side length of each cubic box is L = 8. The FMAP calculations are done on a 80 × 80 × 80 grid, with grid spacing at 0.1. To compensate for discretization errors, the Lennard-Jones potential is adjusted, with ε increased by a factor of 1.06 and σ decreased by a factor of 6/6.06 [9]. (c) Maxwell construction based on the full chemical potential as a function of the particle density. The equal-area chemical potential, indicated by a dotted horizontal line, is calculated using the van der Waals loop in which the excess chemical potential is given by the polynomial fit. The left-most and right-most intersections of the line with the loop are the densities of the dilute and dense phases (indicated by circles).

Fig. 5.

Illustration of Gibbs-ensemble Monte Carlo. Two separated systems are simulated, with the temperature, total particle number, and total volume held constant. (a) Three types of Monte Carlo moves: particle displacement; volume exchange; and particle swap. (b) The densities of the two systems during the simulations at N = 2048 and T = 1.05; the average densities in the second half of the simulations are shown as red horizontal lines. Snapshots of the two systems at the end of 100,000 Monte Carlo (MC) cycles are shown. Each MC cycle consists of, on average, N particle displacements, 5 volume exchanges, and N particle swaps. (c) The binodal, constructed from the average densities of the two systems at different temperatures.

Fig. 6.

Setup and data analysis of slab-geometry simulations. (a) The system, initially in a cubic box, is compressed in all directions to reach a density around 0.7. The smaller cubic box, with side length L, is then padded in the +z and −z directions with empty rectangular boxes that each have length 2L. Finally the system is allowed to equilibrate in molecular dynamics simulations at constant N, V, and T. (b) Calculation of density profile along z. The simulation box (at N = 20000 and T = 0.85) is divided into slabs, with cross section L² and thickness Δz. The number, n_i, of particle in slab i is counted to yield the density ρ_i = n_i/(L²Δz). (c) The density profile (blue curve), averaged over saved snapshots, and its fit to a hyperbolic tangent function (red). (d) Contribution of a pair interaction to the z-dependent pressure tensor, according to Irving and Kirkwood [33]. For every pair of particles (i and j) within the cutoff distance, the number (n_s) of slabs with thickness Δz that span the displacement, z_ij = z_j – z_i, is counted. The contribution of this pair to the pressure tensor (Eq [15]) is then distributed equally among the n_s slabs. Finally, within each slab, the pieces from all the interaction pairs are summed to yield the z-dependent pressure tensor.

Fig. 7.

Comparison of binodals calculated from chemical potentials determined by FMAP and by brute-force insertion. The volumes of the systems at different densities are fixed at L³ with L = 8.

Fig. 8.

Dependence of FMAP binodals on system size. (a) Top: binodals at three system sizes. Relative to the series of systems at L = 8, the particle numbers are doubled and quadrupled, respectively; the volume is increased by nearly the same factors, with L increasing to 10 and 12.8. The grid spacing is fixed at 0.1. The inset shows a zoom into the region near the critical point, where size effect is most prominent. Note the narrowing of the binodals with increasing system size. Bottom: fit of the binodal at L = 12.8 to Eqs [16] and [17]. (b) Effects of system size on the fitting parameters. The critical temperature decreases with increasing L as T_c∞ + C/L² (curve, with T_c∞ = 1.17 and C = 0.9), consistent with the narrowing of the binodals, but the other three parameters lack any trend, as indicated by the comparison against a horizontal line.

Fig. 9.

Details of the three types of Monte Carlo moves. Each Monte Carlo cycle consists of a total of 2N + 5 moves, distributed among the three types with probabilities of N/(2N + 5), 5/(2N + 5), and N/(2N + 5), respectively. (a) For each particle displacement, a particle is randomly selected from the list of N. In the case illustrated, the selected particle is in system I. A displacement is then proposed, with uniform probability within a cube centered at the original position and with side length 2l_max. The would-be change in the total energy, ΔU_I, of system I is calculated, and the proposed displacement is accepted with the probability shown, where ΔU_II = 0. (b) For each volume exchange, the volume of system I is proposed to increase by ΔV, where ΔV is a random value picked between –ΔV_max and ΔV_max with uniform probability, and the volume of system II is proposed to decrease by ΔV. In the case illustrated, ΔV is negative. The proposed volume exchange is accepted with the probability shown. (c) For each particle swap, the system to which a particle is to be inserted is selected between I and II with equal probability. In the case illustrated, system I is selected for particle insertion while system II is chosen for particle removal. The proposed particle swap is accepted with the probability shown.

Fig. 10.

Weak dependence of GEMC binodals on system size. (a) Top: binodals at five system sizes. The average density of the systems of the GEMC simulations is fixed at 0.3, while the total number of particles increases from 128 to 2048. The binodals at the four larger sizes essentially agree with each other, but that at N = 128 shows a minute leftward shift (toward lower densities), especially in the dense-phase branch. (b) Bottom: fit of the binodal at N = 2048 to the law of rectilinear diameters. (b) Effects of system size on the fitting parameters. The values of each parameter at the different system sizes are very close to each other, except for a slightly lower ρ_c and a slightly lower B at N = 128, corresponding to the minute leftward shift and narrowing of the binodal.

Fig. 11.

Illustration of a time step in Langevin dynamics simulations. (a) Motions of particles in one time step. (b) Governing equations of particle motions. (c) Algorithm for advancing particle positions and velocities in one time step. Symbols are: F_i, force on particle i due to interparticle interactions; ζ, damping constant; R_i and R_iα, random force and its Cartesian component; δt, time step; G_i and G_i′, 3-dimensional vectors with components given by random variables with zero mean and unit standard deviation.

Fig. 12.

Dependence of slab-geometry binodals on system size. (a) Top: binodals at six system sizes. Both N and L increase, but the density N/L³ is kept around 0.7. With increasing system size, the binodals broaden near the critical point. Bottom: fit of the binodal at N = 50000 to the law of rectilinear diameters. (b) Dependences of the fitting parameters on system size, modeled by the function q = q_∞ + C/N. For T_c, ρ_c, B, and A, the values of (q_∞, C) are: (1.175, −43.1), (0.316, −16.0), (0.511, 19.8), and (−0.190, −58.5), respectively.

Fig. 13.

Interfacial tension results from slab-geometry simulations at N = 20000. (a) The pressure profile according to the Irving-Kirkwood method (solid curves), and its integration that yields the interfacial tension (dashed curves). (b) Comparison of results by the Kirkwood-Buff method (symbols) and by the Irving-Kirkwood method (lines).

Fig. 14.

Dependence of interfacial tension on system size. (a) Interfacial tension results calculated at four system sizes over a range of temperature. The inset shows a zoom into the data at T = 0.85. (b) Dependence of interfacial tension at T = 0.65 on system size, modeled by the function γ = γ_∞ + D/N^2/3, with (γ_∞, D) = (0.856, 0.801). (c) Corresponding data at T = 0.85, with (γ_∞, D) = (0.468, 0.767). (d) Corresponding data at T = 0.95, with (γ_∞, D) = (0.292, 0.606).

Fig. 15.

Binodals, extrapolated to infinite system size. (a) Extrapolated binoidals (curves), compared with the raw data obtained at the largest system size for each method. The extrapolated binodals are calculated from the law of rectilinear diameters, with the values of T_c, ρ_c, B, and A given by those extrapolated to infinite system size (Figs. 8b, 10b, and 12b). (b) Comparison of the raw data from the three methods at the largest system sizes. (c) Comparison of the binodals extrapolated to infinite system size, by the three methods.

Fig. 16.

Interfacial tension, extrapolated to infinite system size. The filled circles are the γ_∞ values from Fig. 14; the open circle is the value calculated at T = 1.05 and N = 50000.