1. Introduction
This paper presents a systematic study of the asymmetric dumbbell model (ASD) consisting of two different-sized Lennard-Jones (LJ) spheres connected by a rigid bond. The ASD model is designed to be a simple molecular model [
1,
2,
3]. Because of the asymmetry, the liquid phase of the model is easily supercooled, i.e., the model does not readily crystallize [
2]. This makes it suitable for numerical studies as a simple, single-component glass-forming model of non-spherically symmetric constituents [
2,
4,
5,
6,
7,
8].
The focus here is not on supercooling and glass formation, however, but on investigating what happens when the bond length is varied in the less viscous liquid as well as plastic-crystal phases. It was previously shown that for the bond length 0.58 (in units of the largest (A) particle’s radius), the ASD model obeys the hidden-scale invariance symmetry, which manifests itself in strong virial potential-energy correlations in the thermal-equilibrium constant-volume fluctuations [
2]. This implies the existence of so-called isomorphs, which are lines in the thermodynamic phase diagram along which structure and dynamics in reduced units are invariant to a good approximation [
2,
9,
10]. We confirm this below and demonstrate that the existence of isomorphs is robust to bond-length changes. This is done by simulating the ASD model in the liquid phase with bond lengths 0.05, 0.1, 0.2, and 0.5. We also investigate the model’s plastic-crystal phase with a focus on the existence of isomorphs; here the largest bond length simulated is 0.3 (for larger bond lengths the system was liquid at the chosen reference state point).
Our investigation reports data illuminating the degree of isomorph invariance of structure and dynamics. This is done by comparing results along isomorphs with results along isotherms of the same density variation. Specifically, the structure is investigated by monitoring the radial distribution function (RDF), while the dynamics is investigated via the mean-square displacement (MSD) as a function of time, as well as the rotational time-autocorrelation function (RAC). We find good isomorph invariance of both structure and dynamics.
The paper presents the first investigation of the consequences of hidden scale invariance for a system of molecules forming plastic crystals. In the plastic-crystal phase the molecules’ centers are ordered on a crystalline lattice—in the present case a face-centered cubic lattice—while the molecular orientations vary more or less randomly because the molecules are free to rotate [
11,
12,
13,
14,
15]. In this phase we find excellent isomorph invariance, in fact better than in the (isotropic) liquid phase.
2. Model and Simulation Details
The asymmetrical dumbbell (ASD) is a constrained molecular model. It consists of two different-sized Lennard Jones (LJ) spheres, a sphere (A) and a slightly smaller one with lower interaction energies (B). If the two spheres are connected by a rigid bond of length =
, the model mimics toluene [
2,
16]. The values of the parameters used below are as follows. For particle A the distance and energy parameters are
and the mass is
; for particle B one has
,
, and
; for the
AB interaction one has
and
[
16].
We define the (number) density
as the total number of A and B particles (
) divided by the simulation box volume
V;
T is below the temperature. The system is studied by molecular-dynamics simulations in the canonical
ensemble using the Nosé-Hoover thermostat [
17]. The simulated system consists of 4000 ASD molecules in a cubic box with periodic boundaries. Simulations are performed using the open-source Roskilde University molecular dynamics software (RUMD) that runs on GPUs (graphics processing units) [
18]. The constant bond-length constraint was implemented using the method of Ref. [
19].
The leapfrog algorithm is used with time step (in A particle LJ units). First, simulations were carried out for time steps for equilibration and time steps for production runs with time step sizes 0.0005, 0.005, 0.001, 0.002, and 0.01 in order to investigate how constant the energy is at the reference state point defined by and . The bond lengths chosen for these simulations were and . We concluded that the time step in both cases shows a good energy conservation. For the subsequent simulations generating the data, at each state point the simulations again ran for time steps for equilibration and time steps for the production runs.
For the post analysis of structure and dynamics, the units used are the isomorph-theory’s reduced units defined as follows: energy unit: , length unit: , time unit: . Note that these units vary with state point. Exceptions to the use of reduced units apply for density and temperature, which by definition are both unity in reduced units; these quantities, as well as bond lengths, are reported in A particle LJ units.
3. Essential Isomorph Theory
Hidden scale invariance is a symmetry of the potential energy function
in which the configuration
R in terms of the particle coordinates is defined by
. If
and
are two same-density configurations, hidden scale invariance is defined by the following mathematical implication (in which
is a uniform-scaling parameter)
This symmetry is only rigorously obeyed if
is an Euler-homogeneous function plus a constant, but it applies to a quite good approximation for many other systems, e.g., single- and multicomponent LJ systems, the Yukawa and bi-Yukawa pair potentials, the EXP and Morse pair potentials, Mie potentials involving various exponents, etc. [
2,
9,
10,
20,
21,
22]. Equation (
1) may also apply for molecular models like the ASD model. The derivations given below apply in the case of constant-length rigid bonds, i.e.,
not to be scaled with
. Equation (
1) then refers to a uniform scaling of the center-of-mass coordinates with molecular orientations and sizes kept unchanged.
The excess entropy
is defined as the entropy minus that of an ideal gas at the same density and temperature. Any state point of the thermodynamic phase diagram is fully characterized by the two thermodynamic variables
and
T, but other variables may of course also be used to characterize a state point, for instance density and average potential energy. We define the
microscopic excess-entropy function
by [
20]
Here the function
is the excess entropy of the state point
in which
U is the average potential energy. It follows from statistical mechanics that
is proportional to the logarithm of the number of configurations with the same density and potential energy as
R [
23]. Inverting Equation (
2) leads to
in which
is the average potential energy of the state point of density
and excess entropy
.
It can be shown that Equation (
1) implies the function
is scale invariant, i.e.,
[
20]. In this case,
depends only on the configuration’s so-called reduced coordinate vector
(which is of course invariant upon a uniform scaling):
Consequently, Equation (
3) becomes
All identities of the isomorph theory may be derived from Equation (
5) [
20]. In particular, Equation (
5) implies strong correlations between the constant-volume fluctuations of the virial
W and the potential energy
U. Thus, the equilibrium deviations from the average of these two variables,
and
, obey [
24,
25,
26]
Here the so-called density-scaling exponent
is defined and characterized by [
9]
The second (generally valid) identity allows one to calculate from the constant-volume thermal-equilibrium canonical-ensemble fluctuations at the state point in question.
Using the identity
, a first-order Taylor expansion of Equation (
5) leads [
20] to
Consider two state points
and
with average potential energies
and
and the same excess entropy
, and suppose that
and
are equilibrium configurations of the two state points with the same reduced coordinates, i.e.,
. In that case, by elimination of
, Equation (
8) implies that with
and
one has
This implies that
in which
is a constant. Equation (
10) is the original 2009 definition of isomorphic state points [
9], which stated that the canonical probabilities of configurations that scale uniformly into one another are identical along an isomorph (the value of the constant
is irrelevant because probabilities are normalized). This identity of the probabilities of scaled configurations implies that
is constant along an isomorph. In fact, the 2014 version of the isomorph theory that introduced Equation (
1)
defined isomorphs as lines of constant excess entropy in the thermodynamic phase diagram and showed from this that the excess entropy is an isomorph invariant [
20].
It can be shown, either from Equation (
10) [
9] or from Equation (
5) [
20], that the dynamics at two isomorphic state points are identical in the “same movie” sense: Filming the molecules’ motion at one state point results in the same movie at a different, isomorphic state point—except for a uniform scaling of space and time. This implies several dynamic isomorph invariants and that the reduced-unit structure is isomorph invariant in reduced units.
How does one decide whether a given system is expected to have good isomorphs? According to Equation (
6) this is the case when the virial and potential-energy fluctuations are highly correlated, which can be investigated by evaluating their Pearson correlation coefficient,
A system is “R-simple”, i.e., simple in the sense of having good isomorphs and thus an essentially one-dimensional thermodynamic phase diagram, if
at the state points of interest [
20,
24,
25]. This is a somewhat arbitrary criterion, though, in the sense that the property of good isomorph invariance depends on the quantity in question—in practice some reduced-unit quantities are more isomorph invariant than others. Note also that how invariant a given property is, of course, depends on how large a density range that is explored.
4. Results for the Liquid Phase
In this section, we investigate the liquid phase for the bond lengths 0.05, 0.1, 0.2, and 0.5. First, however,
Figure 1a shows a snapshot of the ASD liquid at the bond length previously studied by Ingebrigtsen et al. [
2] (
) [
2] at a liquid state point. The A particles are red, the B particles are blue. We see the typical disorder of a liquid configuration. Reducing the bond length results in a system that is increasingly like the atomic LJ liquid, but even for the smallest bond length (
) we still find typical molecular-liquid behavior (see below).
Isomorphs, i.e., curves of constant excess entropy
, are traced out by numerically integrating Equation (
7). In this way one avoids thermodynamic integration to determine
throughout the phase diagram. Specifically, the last term of Equation (
7) specifies how the density-scaling exponent
may be calculated from an
simulation at the state point in question. Integrating the first-order differential equation defined by the second equality sign in Equation (
7) in order to determine how temperature varies with density along an isomorph, is in principle straightforward. The highly accurate fourth-order Runge–Kutta (RK4) integration method was recently implemented and used for this [
27,
28]. We used that method here with a density change of 1%. For all four bond lengths we started the integration at the reference state point
. Thus, all four isomorphs have this state point in common.
Table 1 gives the pressure, density-scaling exponent, and virial potential-energy correlation coefficient for the four bond lengths at the reference state point.
The isomorphs are not very different (
Figure 1b), but we note that larger bond lengths result in slightly larger temperatures at a given density. This reflects larger density-scaling exponents
for the larger bond lengths, as is evident from
Figure 1c which shows
at the reference state point for several bond lengths. A value of
between 5 and 6 is typical for atomic LJ systems [
25,
26]. In fact, in the liquid phase a density-scaling exponent between 5 and 6 is not only typical of atomic LJ system, but also for diatomic LJ systems that often have values somewhat higher than the atomic case [
29,
30]. Interestingly, density-scaling exponent relatively close to 6 are often found in experiments [
31,
32,
33]. We note a slight drop of
at the largest bond length studied, but the overall
variation is just 15% while, as we shall see, the physics varies considerably when the bond length is changed from 0.05 to 0.5.
As mentioned, the isomorph theory only applies for systems with strong virial potential-energy correlations. To ensure that this requirement is obeyed, we evaluated the correlation coefficient
R (Equation (
11)) along the isomorphs and along
isotherms. The results are shown in
Figure 2 as functions of the density in which the traditional
pragmatic limit for isomorph theory to apply is marked by horizontal dashed lines. We see that the majority of state points are above this line, suggesting that the isomorph theory applies. The strongest correlations are found for the largest bond length studied (0.5).
To investigate to which degree there is isomorph invariance of the reduced-unit structure we performed simulations of the three radial distribution functions (RDF)
defined by the AA, AB, and BB particle pairs [
17] along the isomorphs. The results are shown in
Figure 3, which for comparison also shows data for simulations carried out at the reference-state point isotherm with the same (20%) density variation, i.e., for densities between that of the reference state point, 1.5, and 1.8. All three RDFs are nicely invariant along the isomorphs compared to their variation along the isotherms. We note a slight variance at the first peak maximum. This deviation from isomorph theory has often been observed; it is due to the fact that in the direction of increasing
along an isomorph, the system becomes more hard-sphere like. Consequently, the left hand side of the first peak diminishes when moving along an isomorph in the direction of increasing
. To “compensate” for this effect and maintain a constant overall coordination number (defined as the RDF integral over the first peak), the peak maximum increases. The variance of the first peak maximum has recently been put into a consistent framework by proofs that the so-called bridge function is isomorph invariant [
28].
Focusing henceforth on the isomorph RDFs, we note that the AA RDF does not change very much as the bond length is increased. In fact, this RDF looks much like that of the standard single-component LJ liquid, which is also very similar to that of, e.g., the hard-sphere, Yukawa, and inverse-power-law liquids [
34]. The observed variation of the AB and BB RDFs can be interpreted as a consequence of the near invariance of the AA particle RDF as the bond length is increased: The structure of the ASD liquid may be thought of as primarily determined by the A particles, which behave much like standard LJ particles that are not really bothered by the B particles. This is because the B particles are smaller and have significantly lower interaction energy parameters. In this picture, the B particles are to a significant degree “slaves” of the A particles, and as a consequence one expects what is observed in
Figure 3: The BB RDF changes significantly with increasing bond length, and at bond length 0.5 it is almost constant. This is because the B particles place themselves in many possible positions around the A particles that, as mentioned, more or less have the RDF of a single-component LJ liquid. This effect is present at all bond lengths, but for the smaller bond lengths it does not give rise to any significant spread of the BB RDF because a given B particle is constrained to be close to the A particle of the same molecule. Indeed, the BB particle RDF is very similar to the AA RDF for the small bond length 0.05. This “B slaving A” picture is confirmed by the AB RDF, which (except for the vertical line coming from the intramolecular AB bond) gives data for the AB correlations between different ASD molecules. The AB RDF is also significantly smeared out as the bond length increases, but less so than the BB RDF because the relative order of the A particles is partly inherited by the AB RDF.
Figure 4 shows data for the (reduced) time dependence of the reduced mean-square displacement (MSD) [
17] of the A and B particles in a log-log plot. The data are isomorph invariant to a good approximation, but not invariant along the isotherm except in the short-time ballistic regime where the invariance is a rigorous consequence of the use of reduced units, which implies unity thermal velocity. Focusing on the behavior along the isomorphs, we note that the A particle motion is pretty similar for all bond lengths. This is consistent with the above picture according to which the A particles to a significant extent behave as if the B particles were not present. In the ballistic regime, the B particles move faster than the A particles because of their lower mass. At long times, the A and B particles are constrained to follow each other, resulting in the same long-time MSD for all bond lengths (which, of course, also applies along the isotherms). An interesting feature appears at intermediate times for the B particle MSD, which for small bond lengths have a slight kink that at larger bond lengths develops into a hint of a plateau. We interpret this as an effect of the fact that the fast ballistic motion of the B particles around the A particles eventually “saturates”. This is confirmed in
Figure 5 discussed next, showing that the rotational time-autocorrelation function for short bond lengths has a negative minimum in this range of the reduced times (around 0.1).
Figure 5 shows data for the rotational time-autocorrelation function (RAC) defined as the autocorrelation of the unit vector from particle A to particle B,
, plotted as a function of the reduced time. Because of the normalization, this quantity always starts in unity at time zero. The upper row shows data for the isotherms, the lower for the isomorphs. There is little difference and good collapse in both cases, except for the largest bond length 0.5 for which the data are significantly less invariant along the isotherm. At this bond length there is also a change of behavior in the sense that the RAC never becomes negative as it does at intermediate times for the lower bond lengths. A negative value of the RAC signals a more than 90 degree average rotation of the molecule. This finding for bond length 0.5 is consistent with the above observation that a change of the physics appears to set in at bond lengths 0.3–0.4. It is important to note that, in contrast to the RDF and the MSD, the RAC is not predicted to be isomorph invariant because the moment of inertia is not isomorph invariant (the bond length is fixed and not scaled with the density by the factor
required for a constant reduced-unit moment of inertia).
5. Results for the Plastic-Crystal Phase
We proceed to report results for the plastic-crystal phase for the bond lengths 0.05, 0.1, 0.2, and 0.3 (for bond lengths larger than 0.3 the system was not crystalline at the chosen reference state point of density 2.2 and temperature 0.5).
Table 2 gives the pressure, density-scaling exponent, and virial potential-energy correlation coefficient for the four bond lengths at the reference state point. A typical configuration is shown in
Figure 6. We clearly see that the system is ordered (compare
Figure 1a). A closer inspection reveals that the order is not perfect, however, because the bond directions are disordered. This is a typical signal of a plastic crystal. At a lower temperature, there must be a phase transition to a perfectly ordered phase in which the bond orientations are also ordered, but we did not investigate that transition.
Before proceeding to probe the extent of isomorph invariance, we look at the properties at the reference state point.
Figure 7a shows the AA particle RDFs for the four bond lengths. These show a well-defined order due to the (almost) crystalline ordering of the A particles; the AA particle RDFs are moreover very similar for all bond lengths. In contrast, the AB and BB RDFs vary significantly because of the rotation of the B particles around the A particles. In these cases, the most ordered RDFs are those of the shortest bonds. These observations show that the slaving of the B particles to the A particles (that are well ordered) applies also in the crystalline phase.
Figure 8 shows the MSD of the A and B particles as functions of the reduced time. For the A particles the results are very similar for the different bond lengths, although bond length 0.3 deviates from the three smaller ones by having a larger long-time plateau. We have no good explanation for this, but deviations for the longest bond length from the three others are also noted in some of the later figures. The B particle plateaus vary considerably, with larger plateaus observed for larger bond lengths. This is because the B particles rotate around the A particles and a larger bond length gives them more freedom to do so, resulting in a larger long-time plateau.
Figure 9 shows
and
R at the reference state point as a function of the bond length. We see that
is almost independent of the bond length and has the typical LJ value between 5 and 6 [
25]. The correlation coefficient drops somewhat with increasing bond length, but stays comfortably above 0.9.
As in the liquid phase, we generated isomorphs from the reference state point by integrating Equation (
7) using the RK4 algorithm [
27] for density steps of magnitude 0.01. Additionally, as in the liquid case, we chose to increase the density by 20% from the reference state-point density. This corresponds to realistic density changes of high-pressure experiments. The four isomorphs are shown in
Figure 10. We see that they are quite similar.
Figure 11 shows the reduced RDFs along the isomorphs and, for comparison, along reference-temperature isotherms of the same density variation. For the smallest bond length (0.05), we see a very good collapse of all three RDFs along the isomorph, but not along the isotherm. Moreover, the three RDFs are quite similar, which is a consequence of the fact that all B particles are constrained to be very close to an A particle. For larger bond lengths we see in all cases a very good isomorph collapse and note that, with increasing bond length, the AB and BB particle RDFs become increasingly smeared out compared to the AA particle RDF. The same was seen in the liquid phase (
Figure 3), and as for the liquid phase we interpret this as deriving from rotations of the molecules. In particular, it confirms that the crystals are plastic. Interestingly, the AA particle RDF becomes more invariant along the isotherm as the bond length is increased, but for all bond lengths this quantity is less isotherm than isomorph invariant.
For the MSD as a function of reduced time (
Figure 12), we note that these are all isomorph invariant to a good approximation. In contrast, there is a notable variation along the reference-state-point isotherm for both A and B particles. Only at short times (in the ballistic regime) do we observe invariance along the isotherms, but as mentioned this is a consequence of the definition of reduced units.
Proceeding finally to the RACs,
Figure 13 shows these where the upper two presents the RACs along the isotherms for the four bond lengths and the lower row gives corresponding data for the RACs along the isomorphs. With the notable exception of bond length 0.3, we see a fairly good isomorph invariance of the RACs. Despite the already-mentioned fact that the RAC is not expected to be isomorph invariant because the reduced moment of inertia is not, we find for all bond lengths that the invariance along isomorphs is better than along the isotherms —and much better in the bond-length 0.3 case.