Isomorph invariance of classical crystals’ structure and dynamics

Isomorph invariance of classical crystals’ structure and dynamics

Dan E. Albrechtsen, Andreas E. Olsen, Ulf R. Pedersen, Thomas B. Schrøder, and Jeppe C. Dyre DNRF Center “Glass and Time”, IMFUFA, Dept. of Sciences, Roskilde University, P. O. Box 260, DK-4000 Roskilde, Denmark
August 2, 2019

This paper shows by computer simulations that some crystalline systems have curves in their thermodynamic phase diagrams, so-called isomorphs, along which structure and dynamics in reduced units are invariant to a good approximation. The crystals are studied in a classical-mechanical framework, which is generally a good description except significantly below melting. The existence of isomorphs for crystals is validated by simulations of particles interacting via the Lennard-Jones pair potential arranged into a face-centered cubic (FCC) crystalline structure; the slow vacancy-jump dynamics of a defective FCC crystal is also shown to be isomorph invariant. In contrast, a NaCl crystal model does not exhibit isomorph invariances. Other systems simulated, though in less detail, are the Wahnström binary Lennard-Jones crystal with the Laves crystal structure, monatomic FCC crystals of particles interacting via the Buckingham pair potential and via a novel purely repulsive pair potential diverging at a finite separation, an ortho-terphenyl molecular model, and SPC/E hexagonal ice. Except for NaCl and ice, the crystals simulated all have isomorphs. Based on these findings and previous simulations of liquid models, we conjecture that crystalline solids with isomorphs include most or all formed by atoms or molecules interacting via metallic or van der Waals forces, whereas covalently- or hydrogen-bonded crystals are not expected to have isomorphs. Crystals of ions or dipolar molecules constitute a limiting case for which isomorphs are only expected when the Coulomb interactions are relatively weak. We briefly discuss the consequences of the findings for theories of melting and crystallization.

I Introduction

In many situations the physical properties of crystals are properly described by classical mechanics Ashcroft and Mermin (1976); Kittel (1976); Sidebottom (2012). Of course, the fact that the specific heat vanishes when temperature is lowered towards zero can only be explained by invoking quantization of the phonon field, but otherwise it makes good sense to evaluate a crystal’s structural and dynamic properties from a purely classical description. This is obviously the case for crystals of large particles like those of colloids Bianchi et al. (2012); Juarez et al. (2012) or dust plasmas Lampe et al. (2000); Durniak et al. (2013), but crystals of atoms and molecules not significantly below the melting temperature are also generally well described by classical mechanics Ubbelohde (1965); Bocchetti and Diep (2013). Thus current theories and computer simulations of melting, superheated crystals, etc, are all formulated within a classical, Newtonian framework Forsblom and Grimvall (2005); Bai and Li (2008); Wang et al. (2008); Gallington and Bongiorno (2010); Zhang et al. (2011); Pedersen (2013).

This paper presents computer simulations of crystals of classical-mechanical particles with a focus on the standard Lennard-Jones (LJ) system, demonstrating invariance of structure and dynamics to a very good approximation along the configurational adiabats. Recall that if is the distance between two particles, the LJ pair potential is defined Lennard-Jones (1924) by


This potential diverges as and goes to zero as ; the prefactor ensures that the minimum value of is . For many years the LJ potential has been the standard pair potential for theoretical studies of liquids Andersen (2005); Eaves and Reichman (2009); Hansen and McDonald (2013). It is also the most widely computer-simulated potential, because it has become a standard building block in molecular models and models of complex systems like biomembranes and large biomolecules Mackerell (2004); Stone (2013).

The theory of isomorphs, which was proposed for liquids in Ref. Gnan et al. (2009), predicts that the configurational adiabats in the thermodynamic phase diagram of certain systems are curves of invariance for many properties. The main finding of the present paper is that many solids also have isomorphs; in fact it turns out that the isomorph theory works even better here than for liquids. The theory does not work well for crystalline sodium chloride or for ice, however, which shows that the existence of isomorphs is not a trivial consequence of crystallinity.

Although more simulations are needed, the picture that emerges from the results reported below is that any liquid with isomorphs also has isomorphs in the crystalline phase. Thus “having isomorphs” is a material-specific property that survives a first-order phase transition, an unusual situation in condensed matter physics.

The invariances found along isomorphs have several important consequences. For instance, the fact that the melting curve is itself an isomorph Gnan et al. (2009); Schrøder et al. (2011) implies invariance along this curve of reduced-unit vibrational mean-square displacement, configurational entropy, isochoric specific heat, radial distribution function(s), phonon spectrum, etc. Such melting-line invariances have been known for many years from simulations and experiments Ubbelohde (1965); Born (1939); Tallon (1980); Wallace (2002), but only now get a unified theoretical explanation. In fact, melting line invariants referring to a single phase – liquid or solid – have always presented a challenge to theory because the melting line is where the two phases’ free energies are identical: how can one phase know about the free energy of another? The isomorph theory as validated below for crystalline systems offers a resolution of this paradox for the systems with isomorphs, a class now known as “Roskilde-simple” or just “Roskilde” systems Flenner et al. (2014); Prasad and Chakravarty (2014); Pieprzyk et al. (2014). Further isomorph-theory predictions is that the dynamics of melting, crystal nucleation, and crystal growth in properly reduced units are all invariant along the isomorphs. For a fixed degree of superheating/supercooling this implies temperature/pressure independence of the melting and crystallization rates and mechanisms.

The paper is structured as follows. In Sec. II we give a few details of the computer simulations carried out. Section III presents data for the potential energies of scaled versus original configurations, showing a clear difference in the behavior between the LJ crystal and a NaCl model crystal. For the LJ crystal the results allow one to identify isomorphic state points, and Sec. IV shows data validating the predicted isomorph invariance of structure and dynamics along an isomorph; in contrast, the NaCl crystal does not exhibit these invariances. Section V shows that jumps along an isomorph for the LJ crystal lead to instantaneous thermal equilibration. The findings for the LJ and the NaCl crystals are put into a broader perspective in Sec. VI, which presents data for five other systems, three atomic and two molecular crystals. Section VII gives a discussion.

Ii Some simulation details

All simulations were standard Nose-Hoover simulations carried out using the Roskilde University Molecular Dynamics (RUMD) software package rum (). Periodic boundary conditions were employed throughout. The first system studied below consists of 4000 LJ particles arranged into a face-centered cubic (FCC) lattice rum (). This is the stable crystal structure of the LJ Stillinger (2001) and many other simple systems like, e.g., that defined by the Buckingham potential (Sec. VI). The potential was cut using the shifted-forces method Allen and Tildesley (1987) with cut-off distance independent of the LJ parameter , corresponding to cut-offs larger than at all densities simulated.

For a systems of particles in volume , the density is defined by . The simulations were carried out in the so-called reduced units defined by the length unit , energy unit , and time unit in which is the particle mass. This is useful in order to avoid overflow problems when large density and temperature changes are involved. In practice, the density and temperature were both kept constant equal to unity – instead the pair potential’s energy and length parameters and were adjusted from state point to state point. The same physics is obtained in this way as by varying density and temperature for fixed and because structure and dynamics in reduced units depend only on the dimensionless parameters and . The time step was 0.0025 in reduced units and the time constant of the Nose-Hoover thermostat was kept constant in reduced units.

The NaCl model studied consisted of particles placed in two interpenetrating FCC lattices. The pair potential is a LJ potential plus a Coulomb term Smith and Dang (1994),


Here , , and have different values depending on the three possible atom-atom interactions, see Ref. Smith and Dang, 1994 for details. Due to the long-ranged nature of the Coulomb part of the NaCl potential Kale and Herzfeld (2011), a shifted-force cutoff with was used Hansen et al. (2012).

For the simulations of Sec. VI the crystals were of the following sizes. The Wahnström binary LJ crystal consisted of 864 small and 432 large particles, the Buckingham and the “sum-IPL” crystals were each of 4000 particles, the OTP crystal consisted of 324 molecules, and the ice crystal of 432 molecules.

Iii Potential energies of scaled configurations and the direct isomorph check

This section motivates the isomorph concept by presenting simulation data for how the potential energy of a uniformly scaled crystalline configuration relates to that of the original configuration (taken from an equilibrium simulation). A recent brief review of the isomorph theory and the evidence for it coming from simulations and experiments is given in Ref. Dyre, 2014; the theory of Roskilde systems is given in a series of five comprehensive papers Bailey et al. (2008a, b); Schrøder et al. (2009a); Gnan et al. (2009); Schrøder et al. (2011).

The first step in the investigation of the LJ crystal is to identify the isomorphs in the thermodynamic phase diagram; these are identical to the configurational adiabats Gnan et al. (2009) (see below). Consider a FCC crystal in thermal equilibrium at density and temperature , which in standard LJ units is written as follows: and . The idea is now that for an equilibrium simulation at this state point, each configuration is scaled uniformly to a different density. In Fig. 1(a) the scaling is to the double density, . Each configuration is characterized by the -dimensional configuration vector , and if the original and the scaled configurations are denoted by and , respectively, the two configurations and have the same reduced coordinates and .

Figure 1: (a) and (b) Scatter plots of the potential energies of scaled configurations. (a) Results from simulations of a single-component Lennard-Jones (LJ) face-centered cubic (FCC) crystal at density and temperature . The figure shows a scatter plot of the potential energy per particle when each configuration is scaled uniformly to density , versus the potential energy of the original configuration denoted by at density . Note that upon the compression the potential energies change from small negative to large positive values. The blue dashed line of slope gives the prediction of the repulsive part of the LJ potential (plus a constant); this does not fit the data that have slope 21.14. (b) Similar scatter plot for simulations of a model NaCl crystal Smith and Dang (1994) in which density is also doubled. (c) Plot of virial and potential energy as functions of time in Argon units for the LJ crystal in equilibrium at the state point , demonstrating very strong correlations; (d) similar plot for the NaCl model, showing much weaker correlations. The correlation coefficient given in figures (c) and (d) is defined in Eq. (4).

Figure 1(a) shows a scatter plot of the potential energies of pairs of scaled versus original configuration generated in an equilibrium simulation at . The potential energies of original and scaled configurations correlate strongly. That this is not a trivial effect of the crystal structure and, for instance, a crystal’s approximate harmonic description, is clear from Fig. 1(b) giving similar data for the NaCl crystal model Smith and Dang (1994). The blue dashed line in Fig. 1(a) is the prediction of the repulsive part of the LJ pair potential (displaced vertically by a large amount to fit the data); clearly this term alone cannot explain the strong correlation observed Bailey et al. (2008b). Strong correlations of the potential energies of scaled configurations have been reported for several liquids, including various LJ-type and molecular model liquids Ingebrigtsen et al. (2012a, b), and in fact even for the 10-bead rigid-bond, flexible LJ-chain model Veldhorst et al. (2014).

The linear correlation of Fig. 1(a) implies that a temperature exists such that the following applies. Whenever two configurations and of the densities and , respectively, have the same reduced coordinates, i.e., , one has This implies the isomorph condition Gnan et al. (2009)


The temperature is determined from the best fit line slope of Fig. 1(a), which is , as follows: .

Equation (3) implies (almost) identical canonical-ensemble probabilities of pairs of scaled configurations belonging to the state points and . As a consequence, the configurational (“excess”) entropy of the two state points are (almost) identical Gnan et al. (2009). In fact, Eq. (3) implies that the structure and dynamics of two isomorphic state points are (almost) identical when given in reduced units Gnan et al. (2009).

Two state points are termed isomorphic if they obey Eq. (3) for their physically relevant configurations Gnan et al. (2009). In this way a mathematical equivalence relation is defined in the thermodynamic phase diagram; its equivalence classes are continuous curves termed isomorphs. The above method of identifying isomorphic state points via a scatter plot of scaled potential energies is referred to as the “direct isomorph check” Gnan et al. (2009).

An isomorph is a configurational adiabat because the excess entropy is an isomorph invariant. The opposite does not apply in general, however – all systems have configurational adiabats, but only some systems have isomorphs. The microscopic virial is defined by Allen and Tildesley (1987). It has been shown that a system has isomorphs if and only if it has strong virial potential-energy correlations for its thermal equilibrium constant-volume fluctuations (Appendix A of Ref. Gnan et al., 2009). A pragmatic definition of a Roskilde system is in which is the virial potential-energy (Pearson) correlation coefficient defined Bailey et al. (2008a) by (where is the deviation from the average virial, the same for , and the sharp brackets denote canonical averages)


It follows from Euler’s theorem that systems with a homogeneous potential-energy function, i.e., those for which for some exponent , are the only ones with 100% virial potential-energy correlations and thus exact isomorphs. This means that any realistic system’s isomorph invariants are only approximate. Figures 1(c) and (d) show the normalized virial and potential-energy equilibrium fluctuations as functions of time for the LJ and NaCl crystals. Only the LJ crystal exhibits strong virial potential-energy correlations, so only for this system is the isomorph theory expected to work Gnan et al. (2009).

Figure 2: Radial distribution functions (RDF) of the LJ crystal as functions of the reduced pair distance of Eq. (6) along various lines in the thermodynamic phase diagram. (a, upper subfigure) RDFs calculated for state points along an isomorph, involving more than a factor of two density change. The data collapse demonstrates structural invariance. For comparison, the bottom two figures of (a) give the RDFs from state points along an isochore and an isotherm, respectively, for the same temperature / density variation. (b) A zoom-in on the first peak of the RDF for isomorph (left) and inverse-power-law scaling implying invariance along the line of constant (right), demonstrating that isomorph invariance is not merely a trivial consequences of the repulsive term of the LJ pair potential.

It was recently shown that the properties of Roskilde systems, including the existence of isomorphs, is a consequence of these systems’ “hidden scale invariance” by which is meant the following: two functions of density exist, and , such that for all physically relevant configurations one has Dyre (2013, 2014)


Here is a dimensionless function of , the reduced configuration vector. For a system of particles interacting via the pair potential , the functions and are both given by expressions of the form in which each term corresponds to a term in Dyre (2013).

The physical content of Eq. (5) is that a density change to a good approximation results simply in a linear affine scaling of the potential-energy surface. Since the addition of an overall density-dependent constant to the potential energy does not affect structure and dynamics – though it does, of course, change the free energy and the pressure – the affine scaling of Eq. (5) can be compensated by adjusting the temperature in proportion to , leading to identical canonical probabilities of configurations that scale uniformly into one another Ingebrigtsen et al. (2012c). This is what happens along an isomorph, resulting in reduced-unit invariance of structure and dynamics Gnan et al. (2009); Dyre (2014). Thus Eq. (5) implies the existence of isomorphs. This equation can also be shown to imply the following thermodynamic separation identity in which is the excess entropy per particle Ingebrigtsen et al. (2012c). Since isomorphs are configurational adiabats, this identity implies that the isomorphs are given by the equation We finally note Nagayama (2011); Ingebrigtsen et al. (2012c) that the separation identity is mathematically equivalent to the configuration-space version of the noted Grüneisen equation of state, according to which pressure is a linear function of energy with constants that only depend on density. This is the standard equation of state for solids under high pressure Nagayama (2011).

Figure 3: RDFs of the NaCl crystal for prospective isomorphic state points identified by the direct isomorph check method (Fig. 1(b)). The figure shows the sodium-chloride, chloride-chloride, and sodium-sodium partial RDFs. This crystal does not have isomorphs (at the highest-density state point simulated the systems in fact has melted).

Iv Isomorph invariance of structure and dynamics: Comparing the LJ and sodium chloride crystals

This section presents simulations of some properties predicted to be invariant along an isomorph. The purpose is to validate the existence of isomorphs for the LJ crystal. The LJ simulations are contrasted to simulations of the NaCl crystal for which, based on Fig. 1, one does not expect isomorphs. Prospective isomorphic state points were identified using the direct isomorph-check method described above for simulations at the following reference state points: for the LJ crystal and for the NaCl crystal.

Figure 2(a) shows the radial distribution function (RDF) of the LJ crystal at different state points as a function of the reduced pair distance defined by


The upper subfigure shows RDFs of the LJ crystal along an isomorph, predicted to be invariant. We see that this applies to very good approximation for more than a factor of two density change along the isomorph, which is represented by seven state points. For comparison, the two lower subfigures give the reduced-unit RDFs for the same temperature/density variations keeping density and temperature constant, respectively.

Most of the state points studied correspond to highly compressed solids, so one might guess that the observed structural invariance derives trivially from the repulsive term dominating the LJ potential. If this were the case, however, state points related by should have the same reduced-unit RDF Schrøder et al. (2009a). Figure 2(b) shows a blow up of the first peak of the reduced-unit RDFs along the curve defined by for the same density variation (right subfigure, same reference state point). There is poor collapse compared to the left subfigure, which is a blow up of the Fig. 2(a) isomorph data. Thus the isomorph collapse is not a trivial consequence of the LJ potential’s repulsive term, which was clear already from Fig. 1(a) in which the scatter plot does not have the slope 16 predicted from the term of the LJ potential.

Isomorph invariance of structure is not a general property of crystals. This is evident from Fig. 3, which shows RDFs for prospective isomorphic state points of the NaCl crystal determined by the direct-isomorph-check procedure. No structural invariance is observed. Thus isomorph invariance is not a trivial consequence of the harmonic approximation that generally describes crystals well.

Figure 4: (a) Normalized velocity autocorrelation functions of the LJ crystal along an isomorph, an isochore, and an isotherm (same state points as in Fig. 2). (b) Normalized velocity autocorrelation functions for the NaCl crystal along a prospective isomorph, showing no data collapse (same state points as in Fig. 3).

Since the potential-energy surfaces of isomorphic state points are identical except for a linear, affine scaling (compare the hidden-scale-invariance identity Eq. (5)), not just the structure, but also the dynamics is predicted to be invariant along an isomorph when given in reduced units Gnan et al. (2009); Dyre (2014). For the LJ crystal we checked this by calculating the single-particle velocity autocorrelation function. Figure 4(a) shows the results for the state points of Fig. 2. Good collapse is found for the isomorphic state points (upper subfigure), although the collapse is not quite as good as for the RDFs. The NaCl crystal shows no data collapse (Fig. 4(b)).

Our simulations show that the phonon dynamics of the LJ crystal is isomorph invariant. Perfect crystals like the one simulated have no slow dynamics. In order to test the generality of the isomorph invariance of the dynamics, we investigated also the much slower dynamics of atoms jumping in a defective LJ crystal obtained by removing some particles. The vacancies thus introduced occasionally jump to new lattice positions, a process making atomic diffusion possible that, because it is thermally activated, becomes slow at low temperatures and/or high densities.

Eight particles were removed from a LJ crystal, thus introducing eight vacancies. The effect of vacancy jumps was monitored by evaluating the mean-square displacement (MSD) of the defective crystal’s particles as a function time. The results are shown in Fig. 5 for state points along an isomorph, an isochore, and an isotherm. The isomorph was generated by the direct isomorph check method applied to the defective crystal starting at , leading to an isomorph that is marginally different from one of the perfect LJ crystal. Figure 5 demonstrates that the regime of slow atomic jump dynamics is also isomorph invariant to a good approximation. This is nontrivial, even in view of the isomorph invariance of structure and fast dynamics of perfect crystals demonstrated in Figs. 2 and 4, because the excitation states of vacancy jumps – the barriers to be overcome – have a small canonical probability and contribute little to the direct-isomorph-check plots used to identify isomorphic state points (Fig. 1).

Figure 5: Reduced-unit mean-square displacement (MSD) of the particles of an FCC LJ crystal from which eight particles were removed, making vacancy-jump dynamics possible. The top figure shows the MSD as a function of reduced time for state points along an isomorph identified by the direct-isomorph-check method. The two bottom figures show MSDs calculated for state points along the isochore and isotherm with the same temperature/density variation. Vacancy jump dynamics depends strongly on temperature and density, but along an isomorph these two effects compensate each other to a good approximation.

V Isomorph jumps of the LJ crystal

Because the canonical-ensemble probabilities of scaled configurations of isomorphic state points are identical, a jump between isomorphic state points is predicted to bring the system instantaneously to equilibrium Gnan et al. (2009). Instantaneous equilibration after isomorph jumps has been demonstrated in simulations of supercooled, highly viscous liquids Gnan et al. (2009); Ingebrigtsen et al. (2012b), as well as for the flexible LJ-chain model Veldhorst et al. (2014). Below we test whether this applies also for “isomorph jumps” of perfect LJ crystals, for which the dynamics takes place on the much faster phonon time scale of order picoseconds (in Argon units).

The procedure used for studying a jump in thermodynamic phase space is the following. First one equilibrates the system at one state point. Then one changes the density by scaling all coordinates uniformly, scaling all velocities to the new temperature, and adjusting the temperature of the thermostat to the new value. Finally, one observes whether or not the system relaxes at the new state point by monitoring the time development of a suitable quantity, in casu the potential energy. Figure 6 shows the potential energy per particle after such jumps from three different state points to the same state point. The initial state points were selected to give an isochoric (red), an isothermal (blue), and an isomorphic (green) jump to the final state point. Only the latter shows instantaneous equilibration, i.e., no change of the potential energy after the jump.

Figure 6: Potential energy per particle after jumps from three different state points of the LJ crystal to the state point . Only the jump from a state point isomorphic to the final state point (green curve) leads to instantaneous equilibration. Time is given in Argon units.

Vi Results for five other model crystals

A number of questions arise from our findings for the LJ and NaCl crystals. To answer some of these we present below simulation results for more systems, three atomic and two molecular crystals.

One question is whether complex crystals may also have isomorphs or this property is limited to simple close-packed crystal structures like the FCC lattice. Liquid-state simulations show that the property of strong virial potential-energy correlations – implying isomorphs Gnan et al. (2009) – is “local” in the sense that for Roskilde systems all interactions beyond the first coordination shell are unimportant and may be ignored Ingebrigtsen et al. (2012a). Based on this one would not expect the crystal structure to be important. In order to illuminate this issue we made use of the recently discovered fact Pedersen et al. (2010) that the Wahnström binary LJ system Wahnström (1991) crystallizes into the intricate Laves phase, in which the unit cell contains no less than twelve atoms. Laves phases are observed in some binary metal systems de Graef and McHenry (2007) and, e.g., binary hard-sphere mixtures of certain size ratios crystallize into Laves phases. Laves phases of metallic elements have a number of intriguing properties, for instance they are not plastically deformable at room temperature wik ().

The Laves phase of the Wahnström binary LJ model is characterized as follows Pedersen et al. (2010). The smaller LJ particles are placed in one of two different distorted icosahedral polyhedra, both made up of six small and six large LJ particles. The large particles sit in a 16-vertex coordination polyhedron comprised of twelve small and four large particles. The latter form a hexagonal diamond network where each large-particle neighbor pair shares six small particle neighbors.

Figure 7: (a) Direct isomorph check for the Wahnström binary LJ crystal, which is a Laves phase structure with a unit cell of twelve atoms. (b) All-particle RDF for the two isomorphic state points identified on the basis of the direct isomorph check in (a). The inset shows the crystal structure Pedersen et al. (2010).

Figure 7 shows results from simulations of this system. For details of the LJ parameters we refer to Ref. Wahnström, 1991; our only modification was to consider a perfect crystal of 2/3 small and 1/3 large particles whereas the original Wahnström paper considered a 50-50 composition in the supercooled liquid state. Figure 7(a) shows direct-isomorph-check results for a 36% density increase. Figure 7(b) shows the RDF for all particles (small and large) at two isomorphic state points with the temperature at the high-density state point calculated from the direct isomorph check of Fig. 7(a) (). We see that even complex crystal structures can have isomorphs. We moreover conclude that the reason the NaCl crystal does not have isomorphs is not simply that it is a two-component system.

Figure 8: (a) Direct isomorph check for the Buckingham potential FCC crystal. This system does not have a repulsive IPL term, but instead an exponentially repulsive term (Eq. (7)). (b) RDFs along an isomorph and the corresponding isochore and isotherm, demonstrating isomorph invariance of the structure, much as for the LJ crystals.

Consider next whether crystalline isomorphs occur only for crystals of particles interacting via pair potentials consisting of IPL terms like the LJ potential. To investigate this we studied a FCC crystal composed of particles interacting via the Buckingham pair potential Buckingham (1938),


In contrast to the LJ case this potential’s repulsive term is finite at (leading to a thermodynamic instability because the attractive term diverges at , though for large values of this is no problem in practice). The system simulated is that of Veldhorst et al. (2012). In the unit system defined by and , Fig. 8(a) shows the direct isomorph check for a doubling of the density starting from and from which we in the usual way calculated the temperature of the isomorphic state point at density (). Figure 8(b) shows the RDFs for these two state points supplemented by a third isomorphic state point (upper subfigure). We conclude that for a crystal to have isomorphs it is not necessary that the pair potential is composed of IPL terms.

Figure 9: (a) Direct isomorph check for an FCC crystal of particles interacting via the purely repulsive “sum-IPL” pair potential of Eq. (8). (b) RDFs along an isomorph, isochore, and isotherm, demonstrating isomorph invariance of the structure much as for the LJ systems. (c) The viral potential-energy correlation coefficient (upper subfigure) and the so-called density-scaling exponent Gnan et al. (2009) (lower subfigure) as functions of density along the isomorph simulated in (b). The density-scaling exponent is quite different from two, its value for an IPL potential; this shows that the logarithmic term is important in the density range studied.

A pure IPL system has a homogeneous potential-energy function, 100% virial potential-energy correlations, and perfect isomorphs Gnan et al. (2009). Given that the LJ crystal has isomorphs and that the NaCl crystal does not, one may speculate that the latter system’s “problem” is that it involves not just two IPL terms, but three. In order to test whether the number of IPL terms is crucial for how well isomorph invariance applies, we simulated the novel pair potential that adds infinitely many IPL terms. Carrying out the integration for , the parameter chosen for our simulations, leads to the following purely repulsive “sum-IPL” pair potential


This potential diverges at and can only be studied at densities with average nearest-neighbor distance larger than . We simulated FCC crystals of densities ranging from to in the unit system defined by . The results are reported in Fig. 9, where (a) shows the direct isomorph check for a density change from to . Figure 9(b) shows the RDFs along the isomorph generated from this and two other direct isomorph checks also generated from simulations, covering a factor of three density variation (upper subfigure). The RDFs of the corresponding isochores and isotherms are shown in the lower subfigures.

The sum-IPL system has excellent isomorphs. At low densities the logarithmic term is almost constant and the sum-IPL potential becomes dominated by the term, i.e., approaches a single-IPL potential (that trivially has perfect isomorphs). To show that the IPL term does not dominate at the densities studied, the lower subfigure of Fig. 9(c) gives the density-scaling exponent as a function of density for the four isomorphic state points of the upper subfigure of (b). For the IPL pair potential , so the logarithmic term is clearly important at the densities studied. Thus many power laws are indeed in play here, and we conclude that the presence of several IPL terms does not necessarily imply poor isomorphs; this is not the NaCl system’s “problem”. – The upper subfigure of Fig. 9(c) shows the virial potential-energy correlation coefficient.

Figure 10: RDFs along an isomorph and the corresponding isochore and isotherm of the Lewis-Wahnström OTP crystal in which each molecule consists of three LJ spheres connected by rigid bonds with a angle between the bonds (the RDFs refer to the LJ particles of different molecules). The isomorph was generated in the usual way from direct isomorph checks (not shown); the uniform scaling of the molecules for the direct isomorph checks generating the isomorphs keeps the bond lengths fixed. The figure demonstrates isomorph invariance of the structure, though not as accurate as for the LJ crystal.

The systems studied so far have been atomic crystals. The isomorph theory, however, is not limited to atomic systems. We end the paper by presenting results for two molecular crystals. First, we consider the Lewis-Wahnström ortho-terphenyl (OTP) simple molecular model consisting of three LJ spheres connected by rigid bonds with angle degrees Lewis and Wahnström (1994). This model is difficult to crystallize and easily supercooled. In the crystal the three LJ spheres sit near the sites of a BCC lattice of a crystal with cubatic orientational disorder, i.e., with the molecules aligning randomly along the three Cartesian axis Pedersen et al. (2011). When a system like this in a simulation is scaled to a different density, the rigid bond lengths are kept constant and only the distances between the molecules’ centers of masses are scaled uniformly to the new density.

Figure 10 shows the simulation results for the intermolecular atom-atom RDF along an isomorph of the OTP crystal (upper subfigure) and along the corresponding isochore and isotherm. The isomorph was generated the usual way from direct isomorph checks (not shown). This time the density change is “only” 25% and the collapse is not as good as for the LJ crystal (it is actually similar to that found for the LJ liquid). Nevertheless, this model still exhibits good data collapse along the isomorph. In this connection it should be kept in mind that typical high-pressure experiments often involve only a 5-10% density change, so a change of 25% is already quite large.

Figure 11: Oxygen-oxygen RDF along a prospective isomorph generated from direct isomorph checks (not shown) for a 20% density change of SPC/E hexagonal ice. The figure demonstrates that this crystal does not have isomorphs, a finding that is consistent with the fact that water is known to have almost zero virial potential-energy correlation coefficient Bailey et al. (2008a).

In order to put the OTP results in perspective we finally show results for an ice model. Water has near zero virial potential-energy correlations at ambient conditions (a fact related to water’s famous density maximum at C Bailey et al. (2008a)), so ice is not expected to have isomorphs. We simulated the SPC/E water model Berendsen et al. (1987) with the molecules arranged into a hexagonal crystal lattice. A prospective isomorph was generated from a direct isomorph check plot (showing poor correlations). The result for the oxygen-atom RDF shown in Fig. 11 not surprisingly demonstrates that not all molecular crystals have isomorphs.

Vii Discussion

We have studied the LJ crystal’s structure and dynamics in some detail and demonstrated their reduced-unit invariance along isomorphs, as well as instantaneous equilibration following a jump between two isomorphic state points. The reduced-unit vibrational density of states is also an isomorph invariant, since it is determined by the velocity autocorrelation function Dickey and Paskin (1969). These findings validate hidden scale invariance for the crystalline LJ system as expressed in Eq. (5), a property that was previously studied only for liquids. The simulations show that the isomorph theory, in fact, works even better for the crystalline than for the liquid phase.

It is important to realize that the existence of isomorphs is not just a high-pressure phenomenon – thus the lowest-density state point of the LJ crystal studied ( and ) is close to the triple point ( and ). The isomorph invariances actually continue to negative pressures (results not shown) as long as the system is metastable and does not phase separate into a crystal plus empty space, which happens around density 0.82 depending on the simulation length and system size.

Isomorph invariance does not apply for the NaCl crystal model. Thus the existence of isomorphs in crystals is not a trivial harmonic effect. Why does the NaCl crystal not have isomorphs? One difference between it and the LJ system is that NaCl is a two-component system. Another difference is that the LJ system’s pair potential involves only two IPL terms, whereas the NaCl system’s involves three. Figures 7 and 9 show, however, that neither fact explains the difference, which is probably due to the long-ranged nature of the NaCl system’s strong Coulomb interactions, forces that have been shown to weaken the virial potential-energy correlations in the liquid phase Bailey et al. (2008a); Schrøder et al. (2009b).

For both the LJ and the NaCl systems the crystals and liquids behave in the same way as regards the existence of isomorphs. The same is the case for the five other models studied in Sec. VI. This is not trivial – few non-universal properties survive the first-order transition separating the liquid and solid states of matter. Is it always the case that a liquid with isomorphs solidifies into a crystal that also has isomorphs? This is an open question, but our simulations so far indicate that the answer is probably yes because the crystal has always turned out to have stronger virial potential-energy correlations than the liquid. Of course, this still makes possible the interesting situation of a crystal with isomorphs that melts into a liquid without.

That the existence of isomorphs in crystals is not limited to the single-component LJ system’s simple FCC crystal is illustrated by the other crystalline systems studied briefly in Sec. VI. In view of the simulation results of this paper and previous experiments and liquid-state simulations, in particular on supercooled liquids Ingebrigtsen et al. (2012a), there is good reason to believe that most or all metals and van der Waals bonded crystals have isomorphs and that this also applies for weakly ionic or dipolar crystals. On the other hand, covalently, hydrogen-bonded, and strongly ionic or dipolar crystals are not expected to have isomorphs. More work is needed to get a full overview of the situation, however.

A recent reformulation of the isomorph theory Schrøder and Dyre () throws new light on the condition for a crystal to have isomorphs. In the new formulation a Roskilde system is characterized by the property that the order of the potential energies of configurations at one density is maintained when these are scaled uniformly to a different density. If and are two physically relevant configurations at one density, this translates into the requirement . As shown in Ref. Schrøder and Dyre, this condition implies all the fundamental characteristics of Roskilde systems Gnan et al. (2009); Dyre (2014), including the existence of isomorphs identical to the configurational adiabats, invariance of structure and dynamics along the isomorphs, and strong virial potential-energy correlations for constant-density fluctuations. To relate this new characterization of Roskilde systems to crystals, we assume that the crystal’s potential energy is well described in the harmonic approximation. Adopting furthermore the Debye approximation, there are just two state-point dependent quantities characterizing the phonon spectrum, the transverse and longitudinal sound velocities. If the ratio between these is state-point independent, it is easy to see that the above inequality is obeyed because all phonon modes scale in the same way when density is changed. A more general formulation is the following: if the harmonic approximation applies and all phonon modes have the same Grüneisen parameter , the crystal is a Roskilde system.

It has been conjectured that at sufficiently high pressure all liquids have strong virial potential-energy correlations and thus isomorphs Papini et al. (2011). If this is confirmed, the same presumably applies for crystalline solids, implying that all crystals have isomorphs at sufficiently high pressure. This may explain the seemingly universal applicability of the high-pressure Grüneisen equation of state (expressing a linear relation between pressure and energy with a constants that are a function only of the density) Nagayama (2011), which is equivalent to the thermodynamic separation identity characterizing a system with isomorphs, which in its configuration-space version has been shown to follow from Eq. (5) Nagayama (2011); Ingebrigtsen et al. (2012c).

Crystals with isomorphs have a thermodynamic phase diagram that is for many purposes effectively one- instead of two-dimensional. For a system with isomorphs the melting and freezing lines are themselves both isomorphs Pedersen (2013); Gnan et al. (2009); Schrøder et al. (2011). This is because the Boltzmann probabilities for crystalline configurations are isomorph invariant, implying that an isomorph cannot cross the melting or freezing lines. Note that this is a topological argument with no reference to free energies. As a consequence, crystals of Roskilde systems, i.e., systems with hidden scale invariance and obeying Eq. (5), are predicted to have constant excess entropy and invariant reduced-unit phonon spectra along the melting line. For such solids it is the excess entropy rather than the free energy that governs the reduced-unit structure and dynamics along the melting line. This fact resolves the apparent paradox mentioned in the Introduction of melting-line invariants referring only to a single phase (liquid or solid) Nieves and Sinno (2011); Grimvall et al. (2012); Lemarchand (2013). For instance, the Lindemann melting criterion – according to which a crystal melts when the vibrational mean-square displacement reaches a certain fraction of the interatomic distance – is predicted to be invariant along the melting line isomorph, i.e., melting takes place for the same reduced vibrational mean-square displacement. Such freezing/melting line invariances have been reported for LJ crystals Luo et al. (2005); Chakravarty et al. (2007) and in experiments Ubbelohde (1965), though occasionally with some deviations at the lowest pressures.

There are only few experimental studies of how the nucleation rate and mechanism and how the crystal growth rate and mechanism vary with pressure and temperature for a supercooled liquid Frekel and McTague (1980); Kelton (1991); Auer and Frenkel (2001). The melting and freezing lines are isomorphs in parallel to a series of isomorphs in the coexistence as well as the supercooled liquid phases. Consequently, the isomorph theory predicts that nucleation and crystal growth properties depend only on the degree of supercooling as quantified by the excess entropy, not separately on pressure or temperature. This is predicted to apply for metallic, van der Waals bonded, and weakly ionic or dipolar supercooled liquids, but not for covalently bonded, hydrogen-bonded, or strongly ionic or dipolar supercooled liquids.

This paper has demonstrated that the existence of isomorphs is not limited to the standard, isotropic liquid state. One may ask whether isomorphs are present also in other anisotropic fluid states or, e.g., when a liquid is confined to certain geometries or subjected to external fields. Although more work is needed in this regard, there are strong indications that the answer is in the affirmative. For instance, simulations have shown that isomorphs exist for some liquids under nanoconfinement Ingebrigtsen et al. (2013), as well as for liquids undergoing a linear or non-linear shear deformation described by the so-called SLLOD equations of motion Separdar et al. (2013). In reference to experiments, consider the transition between the nematic and isotropic phases of a liquid crystal. For a system with isomorphs, just as for melting, the nematic-isotropic transition line is an isomorph. Indeed, it has been shown for several liquid crystals that in the two-dimensional thermodynamic phase diagram the molecular “flip-flop” reorientational relaxation time is invariant along the nematic-isotropic phase transition line Roland (2008); Urban and Roland (2011). In other words, while the transition temperature varies with pressure along the transition line, the reorientational time does not. This is what one expects from isomorph invariance since the transition line is an isomorph along which the (reduced) relaxation time is consequently invariant.

At temperatures much below melting the structure, dynamics, and thermodynamics of a crystal become increasingly dominated by quantum effects. An important question for future work is whether there also in this region of the thermodynamic phase diagram are simplifying features for crystals of Roskilde systems, i.e, those that have classical-mechanical isomorphs at higher temperatures.


  • Ashcroft and Mermin (1976) N. W. Ashcroft and N. D. Mermin, Solid State Physics (Holt, Rinehart and Winston, 1976).
  • Kittel (1976) C. Kittel, Introduction to Solid State Physics, 5th Ed. (Wiley, New York, 1976).
  • Sidebottom (2012) D. L. Sidebottom, Fundamentals of Condensed Matter and Crystalline Physics (Cambridge University Press, 2012).
  • Bianchi et al. (2012) E. Bianchi, G. Doppelbauer, L. Filion, M. Dijkstra, and G. Kahl, J. Chem. Phys. 136, 214102 (2012).
  • Juarez et al. (2012) J. J. Juarez, S. E. Feicht, and M. A. Bevan, Soft Matter 8, 94 (2012).
  • Lampe et al. (2000) M. Lampe, G. Joyce, G. Ganguli, and V. Gavrishchaka, Phys. Plasmas 7, 3851 (2000).
  • Durniak et al. (2013) C. Durniak, D. Samsonov, J. F. Ralph, S. Zhdanov, and G. Morfill, Phys. Rev. E 88, 053101 (2013).
  • Ubbelohde (1965) A. R. Ubbelohde, Melting and Crystal Structure (Clarendon, Oxford, 1965).
  • Bocchetti and Diep (2013) V. Bocchetti and H. T. Diep, J. Chem. Phys. 138, 104122 (2013).
  • Forsblom and Grimvall (2005) M. Forsblom and G. Grimvall, Nat. Mater. 4, 388 (2005).
  • Bai and Li (2008) X.-M. Bai and M. Li, Phys. Rev. B 77, 134109 (2008).
  • Wang et al. (2008) N. Wang, S. I. Rokhlin, and D. F. Farson, Nanotechnology 19, 415701 (2008).
  • Gallington and Bongiorno (2010) L. C. Gallington and A. Bongiorno, J. Chem. Phys. 132, 174707 (2010).
  • Zhang et al. (2011) H. Zhang, P. Kalvapalle, and J. F. Douglas, J. Phys. Chem. B 115, 14068 (2011).
  • Pedersen (2013) U. R. Pedersen, J. Chem. Phys. 139, 104102 (2013).
  • Lennard-Jones (1924) J. E. Lennard-Jones, Proc. R. Soc. London A 106, 441 (1924).
  • Andersen (2005) H. C. Andersen, Proc. Natl. Acad. Sci. (USA) 102, 6686 (2005).
  • Eaves and Reichman (2009) J. D. Eaves and D. R. Reichman, Proc. Natl. Acad. Sci. (USA) 106, 15171 (2009).
  • Hansen and McDonald (2013) J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids: With Applications to Soft Matter (Academic, New York, 2013), 4th ed.
  • Mackerell (2004) A. D. Mackerell, J. Comput. Chem. 25, 1584 (2004).
  • Stone (2013) A. Stone, The Theory of Intermolecular Forces (Oxford University Press (Oxford), 2nd Edition, 2013).
  • Gnan et al. (2009) N. Gnan, T. B. Schrøder, U. R. Pedersen, N. P. Bailey, and J. C. Dyre, J. Chem. Phys. 131, 234504 (2009).
  • Schrøder et al. (2011) T. B. Schrøder, N. Gnan, U. R. Pedersen, N. P. Bailey, and J. C. Dyre, J. Chem. Phys. 134, 164505 (2011).
  • Born (1939) M. Born, The Journal of Chemical Physics 7, 591 (1939).
  • Tallon (1980) J. L. Tallon, Phys. Lett. A 76, 139 (1980).
  • Wallace (2002) D. C. Wallace, Statistical Physics of Crystals and Liquids (World Scientific, Singapore, 2002).
  • Flenner et al. (2014) E. Flenner, H. Staley, and G. Szamel, Phys. Rev. Lett. 112, 097801 (2014).
  • Prasad and Chakravarty (2014) S. Prasad and C. Chakravarty, J. Chem. Phys. 140, 164501 (2014).
  • Pieprzyk et al. (2014) S. Pieprzyk, D. M. Heyes, and A. C. Branka, Phys. Rev. E 90, 012106 (2014).
  • (30) All simulations were performed using a molecular dynamics code optimized for NVIDIA graphics cards, which is available as open source code at
  • Stillinger (2001) F. H. Stillinger, J. Chem. Phys. 115, 5208 (2001).
  • Allen and Tildesley (1987) M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Oxford Science Publications, 1987).
  • Smith and Dang (1994) D. E. Smith and L. X. Dang, J. Chem. Phys. 100, 3757 (1994).
  • Kale and Herzfeld (2011) S. Kale and J. Herzfeld, J. Chem. Theory. Comput. 7, 3620 (2011).
  • Hansen et al. (2012) J. S. Hansen, T. B. Schrøder, and J. C. Dyre, J. Phys. Chem. B 116, 5738 (2012).
  • Dyre (2014) J. C. Dyre, J. Phys. Chem. B 114, DOI:10.1021/jp501852b (2014).
  • Bailey et al. (2008a) N. P. Bailey, U. R. Pedersen, N. Gnan, T. B. Schrøder, and J. C. Dyre, J. Chem. Phys. 129, 184507 (2008a).
  • Bailey et al. (2008b) N. P. Bailey, U. R. Pedersen, N. Gnan, T. B. Schrøder, and J. C. Dyre, J. Chem. Phys. 129, 184508 (2008b).
  • Schrøder et al. (2009a) T. B. Schrøder, N. P. Bailey, U. R. Pedersen, N. Gnan, and J. C. Dyre, J. Chem. Phys. 131, 234503 (2009a).
  • Ingebrigtsen et al. (2012a) T. S. Ingebrigtsen, T. B. Schrøder, and J. C. Dyre, Phys. Rev. X 2, 011011 (2012a).
  • Ingebrigtsen et al. (2012b) T. S. Ingebrigtsen, T. B. Schrøder, and J. C. Dyre, J. Phys. Chem. B 116, 1018 (2012b).
  • Veldhorst et al. (2014) A. A. Veldhorst, J. C. Dyre, and T. B. Schrøder, J. Chem. Phys. 141, (in press) (2014).
  • Dyre (2013) J. C. Dyre, Phys. Rev. E 88, 042139 (2013).
  • Ingebrigtsen et al. (2012c) T. S. Ingebrigtsen, L. Bøhling, T. B. Schrøder, and J. C. Dyre, J. Chem. Phys. 136, 061102 (2012c).
  • Nagayama (2011) K. Nagayama, Introduction to the Grüneisen Equation of State and Shock Thermodynamics (Amazon Digital Services, Inc., 2011).
  • Pedersen et al. (2010) U. R. Pedersen, T. B. Schrøder, J. C. Dyre, and P. Harrowell, Phys. Rev. Lett. 104, 105701 (2010).
  • Wahnström (1991) G. Wahnström, Phys. Rev. A 44, 3752 (1991).
  • de Graef and McHenry (2007) M. de Graef and M. E. McHenry, Structure of Materials (Cambridge University Press, 2007).
  • (49)
  • Buckingham (1938) R. A. Buckingham, Proc. R. Soc. London A 168, 264 (1938).
  • Veldhorst et al. (2012) A. A. Veldhorst, L. Bøhling, J. C. Dyre, and T. B. Schrøder, Eur. Phys. J. B 85, 21 (2012).
  • Lewis and Wahnström (1994) L. J. Lewis and G. Wahnström, Phys. Rev. E 50, 3865 (1994).
  • Pedersen et al. (2011) U. R. Pedersen, T. S. Hudson, and P. Harrowell, J. Chem. Phys. 134, 114501 (2011).
  • Berendsen et al. (1987) H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269 (1987).
  • Dickey and Paskin (1969) J. M. Dickey and A. Paskin, Phys. Rev. 188, 1407 (1969).
  • Schrøder et al. (2009b) T. B. Schrøder, U. R. Pedersen, N. P. Bailey, S. Toxvaerd, and J. C. Dyre, Phys. Rev. E 80, 041502 (2009b).
  • (57) T. B. Schrøder and J. C. Dyre, arXiv:1406.2216 (2014).
  • Papini et al. (2011) J. J. Papini, T. B. Schrøder, and J. C. Dyre, arXiv:1103.4954 (2011).
  • Nieves and Sinno (2011) A. M. Nieves and T. Sinno, J. Chem. Phys. 135, 074504 (2011).
  • Grimvall et al. (2012) G. Grimvall, B. Magyari-Koepe, V. Ozolins, and K. A. Persson, Rev. Mod. Phys. 84, 945 (2012).
  • Lemarchand (2013) C. A. Lemarchand, The Journal of Chemical Physics 138, 034506 (2013).
  • Luo et al. (2005) S.-N. Luo, A. Strachan, and D. C. Swift, J. Chem. Phys. 122, 194709 (2005).
  • Chakravarty et al. (2007) C. Chakravarty, P. G. Debenedetti, and F. H. Stillinger, J. Chem. Phys. 126, 204508 (2007).
  • Frekel and McTague (1980) D. Frekel and J. P. McTague, Ann. Rev. Phys. Chem. 31, 491 (1980).
  • Kelton (1991) K. F. Kelton, Solid State Phys. 45, 75 (1991).
  • Auer and Frenkel (2001) S. Auer and D. Frenkel, Nature 409, 1021 (2001).
  • Ingebrigtsen et al. (2013) T. S. Ingebrigtsen, J. R. Errington, T. M. Truskett, and J. C. Dyre, Phys. Rev. Lett. 111, 235901 (2013).
  • Separdar et al. (2013) L. Separdar, N. P. Bailey, T. B. Schrøder, S. Davatolhagh, and J. C. Dyre, J. Chem. Phys. 138, 154505 (2013).
  • Roland (2008) C. M. Roland, Soft Matter 4, 2316 (2008).
  • Urban and Roland (2011) S. Urban and C. M. Roland, J. Non-Cryst. Solids 357, 740 (2011).
Ulf R. Pedersen gratefully acknowledges the support of the Villum Foundation’s grant VKR-023455. The center for viscous liquid dynamics “Glass and Time” is sponsored by the Danish National Research Foundation via grant DNRF61.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description