Reentrant melting and multiple occupancy crystals of bounded potentials: Simple theory and direct observation by molecular dynamics simulations

Aspects of the phase coexistence behavior of the generalized exponential model (GEM-m) and bounded versions of inverse power potentials based on theory and molecular dynamics (MD) simulation data are reported. The GEM-m potential is $\varphi \left(r\right)=exp(-r^m)$, where $r$ is the pair separation and $m$ is an adjustable exponent. A simple analytic formula for the fluid-solid envelope of the Gaussian core model which takes account of the known low- and high-density limiting forms is proposed and shown to represent the simulation data well. The bounded inverse power (BIP) potential is $\varphi \left(r\right)=1/{(a^q+r^q)}^{n/q}$, where $a$, $n$, and $q$ are positive constants. The BIP potential multiple occupancy crystal or cluster crystals are predicted to form when $q>2$ and $a>0$, for $n>3$, which compares with the corresponding GEM-m condition of $m>2$. Reentrant melting should occur for the BIP potential when $q\le 2$ and $a>0$. MD simulations in which the system was gradually compressed at constant temperature using the BIP potential produced cluster states in the parameter domain expected but it was not possible to establish conclusively whether a multiply occupied crystal or a cluster fluid had formed owing to assembly structural fluctuations. The random phase approximation reproduces very well the BIP MD energy per particle without any discontinuities at the phase boundaries. The Lindemann melting rule is shown analytically to give a more rapidly decaying reentrant melting curve boundary than the so-called melting indicator (MI) empirical melting criterion which has also been investigated in this study. The MI model gives a better match to the high-density phase boundary for small $m$ and $q$ values.

Figure 1

(a) The general exponential potential (GEM-m) defined in Eq. (1) for several $m$ values which are indicated in the figure. (b) The BIP potential given in Eq. (2). In the examples shown the bounding parameter, $a=1$, and the values of $n$ and $q$ are given in that order in the figure. Data for $n=18$ are given in the left-hand panel (shifted up by 1) and the $n=6$ case is presented in the right-hand panel. The thick black lines on both frames are for the GCM potential which acts a typical reference bounded potential.

Figure 2

(a) The fluid-solid phase diagram line of the GCM taken from Ref. [17] (“MD Pr.”), and Mausbach et al. [55] (“MS”). The green curve (MD Pr.) is the bcc/fcc boundary line from Ref. [17]. The system is a fluid at all temperatures above about 0.009. Direct phase change MD simulations from this work are denoted by “MD TW,” and “MI” is Eq. (15) with $L=190$. (b) The phase diagram of the exponential potential. The magenta line is the LMR prediction. The solid blue squares are results from MD simulations [18]. The lines are the Lindemann (bcc) LMR formula and the melting indicator (MI) defined in Eq. (15) with $L=190$. The fit formula in Eq. (14) denoted by “CM” is hardly distinguishable from the MI curve with these parameters.

Figure 3

(a) The dependence of $T_{\max }$ on $m$ for the GEM-m potential class. Data calculated numerically at intervals of 0.25 for $m$ using the LMR and CM formulas. (b) The same as in (a) except that the dependence of ${\rho }_{\max }$ on $m$ is given. The polynomial fit parameters for this quantity are given in Table 2.

Figure 4

(a) A comparison of the $T\left(\rho \right)$ fluid-solid boundary lines from Eqs. (9) and (10) for LMR and Eq. (15) for MI. The parameters used were for a bcc lattice with $L=200$, for the GEM-m potential with $m=1$ and $m=2$. The figure confirms that the high-density boundary passes through $T=0$ at a lower density for the LMR formula, as discussed in the main text in regard to Eq. (17). (b) The same as for (a) except using the curves for the BIP potential, i.e., where Eq. (21) is used in Eq. (9) (LMR) and Eq. (15) was taken for the MI method.

Figure 5

The radial distribution functions of the exponential potential, $exp(-r^m)$, where $m=2$ (i.e., GCM) and $T=0.002$. (a) The RDFs for densities increasing from bottom to top between $\rho =0.12$ and 1.2 in 50 equal density intervals. The reentrant melting density is about $0.67\pm 0.02$. (b) The same as for (a) except the density is decreased in steps in reverse order. Again the lowest-density RDF is shown at the bottom. The system had 3456 particles and was started at low density from a bcc lattice.

Figure 6

As for Fig. 5 except that $m=3$, and the starting lattice was a bcc crystal. The temperature is 0.004. The density range is between 0.05 and 1.2. The system transforms from a bcc crystal to a multiple occupancy state at a density of $0.71\pm 0.02$ with $N=3456$.

Figure 7

The estimated phase diagrams for (a) $m=2.5$ and (b) $m=3.0$. The vertical red line and the solid blue circles indicate the MOP transition obtained from the simulations of this work. The sloping linear magenta line on the left-hand side of the figure denotes the fluid-solid MOC transition temperature from the MFA approximation from Eq. (22). For $m=2.5$ and 3, $\tilde{\varphi }\left(k_{*}\right)$ is $-0.02908$ and $-0.06672$, respectively, employed in Eq. (22). The corresponding $k_{*}$ are 5.60 and 5.33, respectively. The black horizontal line indicates the maximum in the fluid-solid envelope of this work obtained directly by MD.

Figure 8

(a) The density dependence of the fluid-solid coexistence temperature, $T$, of the inverse power [i.e., the BIP potential from Eq. (2) with $a=0$]. The values of $n$ are given in the figure. Continuous curves using IP scaling formulas benchmarked with Monte Carlo simulation data from Refs. [30, 31] are shown. These are $T={\rho }^{4/3}/10.2070$, $T={\rho }^{6/3}/5.4636$, $T={\rho }^{12/3}/1.9659$, and $T={\rho }^{18/3}/1.2929$ for $n=4$, 6, 12, and 18, respectively. Also given as symbols are the values from the Lindemann formulas (LMR) in Eqs. (9) and (10). The bcc lattice formula gave the best agreement with the simulation-derived equation of state data (EoS) using LMR on taking $r_{nn}={(\alpha 3^{3/2}/4\rho )}^{1/3}$ where the factor $\alpha =0.95$ to improve agreement with the EoS. (b) The same as for (a) except the corresponding BIP potential curves are presented adopting variable $a$ (given in the figure) using $n=6$ and $q=2$ in each case.

Figure 9

As for Fig. 8 except (a) the LMR coexistence curve for the BIP potential systems with $n=6$ and $q=4$ with variable $a$ are shown. (b) The same as for (a) except $q=2$, $a=1$, and the indicated BIP exponent, $n$, is varied.

Figure 10

(a) Radial distribution functions plotted as a function of $r{\rho }^{1/3}$ for the BIP potential in the density range 0.5–20.0. Density increases from bottom to top in the figure. The simulations started from a bcc lattice and $N=3456$. The parameters used were $a=1$, $n=12$, and $q=2$. (b) The variation in the values of certain properties as a function of density during the increasing density part of the MD cycle for the simulation of (a). The height of the first peak in $g\left(r\right)$, or “gr-max,” is shown. The ratio of the first RDF minimum to the first maximum, or “gr-min/max,” is also given. The average potential energy per particle, $u$, is presented as open circles. The solid blue line coincident with the $u$ data is from the RPA formula in Eq. (30), where $\tilde{\varphi }\left(0\right)=0.539744$.

Figure 11

As for Fig. 10, except $q=2.1$ and $\tilde{\varphi }\left(0\right)=0.620338$ which was obtained by numerical integration. The other parameters used were $T=0.002$, $N=3456$, $a=1$, and $n=12$. (b) The average occupancy number, $n_c$, is also presented as $n_c-1$. The particle number density range covered was from 0.5 to 25.

Figure 12

As for Fig. 11, except $q=2.2$ and the position of the first peak in $g\left(r\right)$, or “rmax,” is shown in (b). The parameters used were $T=0.002$, $N=3456$, $a=1$, and $n=12$. The density range is from 0.5 to 50. The RPA constant for Eq. (30) is $\tilde{\varphi }\left(0\right)=0.703419$.

Figure 13

As for Fig. 12, except $q=2.3$ and the position of the first peak in $g\left(r\right)$, or “rmax,” is shown in (b). The density range is from 0.5 to 50, and the other parameters are $N=6750$, $T=0.002$, $a=1$, and $n=12$. The RPA constant for Eq. (30) is $\tilde{\varphi }\left(0\right)=0.788368$.

Figure 14

Perspective snapshot particle assembly images showing clustering. (a) GEM with $T=0.4$, $\rho =7.0$, $m=4.0$, and $N=3456$. (b) BIP with $T=0.002$, $n=12$, $q=2.20$, $a=1.0$, and $N=3456$.