I. Introduction

The strong, complicated, and continual interactions that operate in condensed phases determine both static and dynamic properties of those material systems. Establishing quantitative, predictive connections between knowledge of those interactions on one hand, and the observable macroscopic properties on the other hand, remains a major challenge to condensed matter theory. In particular, this is true for the self-diffusion process, one of the most basic attributes of molecular motion. This paper has as its primary goal the clarification of diffusive motion from the viewpoint of the multidimensional potential energy "landscape" representation,^1,2 whereby steepest descent basins, and their embedded "inherent structures", serve as descriptive analytical tools.^3-5 To illustrate our theoretical ideas, we have carried out, and report results from, a series of molecular dynamics computer simulations. These simulations utilize the classical single-component Dzugutov model,^6-8 over a wide temperature range. We are pleased to note that Keyes and Chowdhary have carried out an analogous investigation for the distinctly different Lennard-Jones model,⁹ and any similarities and contrasts between the results obtained for these two cases could provide a basis for enhanced understanding of molecular details of diffusive processes. The introductory section to ref 9 contains a detailed general discussion of the physical insights provided by the potential landscape/inherent-structure approach.

The following section briefly reviews the basic formulas describing diffusive motion, the procedure for exhaustive division of the multidimensional potential energy landscape into basins, and the use of inherent structures (potential energy minima) as a descriptive approach. Section III reviews the Dzugutov model, explains why it was chosen for the present study, and describes how it has been implemented in our molecular dynamics simulation. Section IV presents numerical results that have emerged from those simulations. Conclusions and discussion in section V end our presentation.

II. Background Formalism

Throughout this paper, attention will focus on classical N-body systems occupying volume V, in a state of equilibrium that is specified by the intensive variables temperature T and particle number density = N/V. Diffusive particle motions that are driven by local dynamics can be represented by the mean-square displacement, versus elapsed time t, for any given "tagged" particle j (1 j N), averaged over the equilibrium ensemble. In the long-time asymptotic limit this mean-square displacement rises at a rate proportional to the self-diffusion constant:¹⁰

Although a best linear fit to F(t) in its asymptotic regime has a slope specified by eq II.1, that linear fit does not generally extrapolate back to the origin. Instead, it will typically exhibit a vertical offset:

An alternative, but closely related, representation of the diffusive process emerges naturally from the steepest descent mapping operation on the multidimensional potential energy hypersurface.^2-5 This operation connects any N-particle configuration r₁...r_N (with zero-measure exceptions) to that of a nearby minimum, to be denoted by R₁...R_N. The latter is conventionally called an "inherent structure". The set of all configurations that map onto the same inherent structure defines a steepest descent basin for that inherent structure. The set of all basins tiles (i.e., exhaustively covers) the entire configuration space. Although each basin contains only a single minimum by construction, it can in principle contain one or more embedded saddles of various orders. Basins are topologically connected but need not be singly connected.

The Newtonian trajectory r_j(t) for any particle j can be resolved into its inherent-structure component R_j(t), plus an intrabasin displacement S_j(t):

Accumulated experience has shown that interbasin transitions, viewed from the respective inherent-structure (IS) configurations, typically entail localized rearrangements of essentially onlyO(1) particles out of the full set of N, and furthermore those localized displacements involved are themselves only O(1) in magnitude.^11-13 This implies that F_IS(t), the inherent-structure mean-square-displacement analogue of F(t),

At low temperature, where the equilibrium phase of the system will be a crystal, or where the system may have been quenched into a low temperature glass, dynamical motions will predominantly be small amplitude harmonic vibrations about a single inherent structure. In this circumstance F_IS(t) will remain equal to zero for any reasonable observation time, but F(t) will rise quickly above zero due to those intrabasin vibrations before leveling out to the value C(T,) > 0. At somewhat higher, but still low, temperatures both F_IS and F will rise slowly with increasing t on account of the presence of slow diffusion, but it is reasonable to anticipate that at intermediate to long times

On account of the piecewise-constant character of the function R_j(t), the short-time behavior of F_IS(t) differs qualitatively from that of F(t). Specifically, the leading order in t is linear:

Equations II.10 and II.11 are an appropriate description for systems with modest, finite numbers of particles such as that chosen for our simulation in sections III and IV. However, it has to be realized that in considering very large systems, and the thermodynamic large-system limit, the interbasin transition rate is an extensive quantity.^11,14 This stems from the fact that the number of locations at which a localized particle rearrangement could occur at any instant is asymptotically proportional to the system size. Therefore, in a very large system the most likely first transition after t = 0 will be remote from the tagged particle j, and will have only an extremely small effect on j due to medium elasticity. Consequently, the t 0 behavior of the infinite-system-limit F_IS(t) may not be simply a linear form indicated in eq II.10 for the finite-N version of that function. This limiting behavior is a subtlety whose detailed analysis lies beyond the scope of the present paper.

It should not be overlooked that comparisons between Newtonian positions r_j(t) and the inherent-structure positions R_j(t) permit evaluation of the instantaneous root-mean-square return distance

III. Simulation Procedure

To illustrate concepts outlined in previous section II with computer simulation, it is advantageous to select a classical model that is relatively simple, that possesses a substantial liquid-state temperature range, and that affords a useful contrast to the Lennard-Jones model that has been previously considered.⁹ With these criteria in mind, we have chosen the Dzugutov pair-potential model for study. When it is expressed in the usual reduced energy and length units, this pair potential has the following form:⁶

Our liquid-phase numerical simulations involved N = 256 particles, confined to a cubic unit cell with a volume chosen to fix the reduced density at = 0.85. Periodic boundary conditions applied. Particle mass was set equal to unity, and the energy and length units used in eqs III.1-III.4 above served to define reduced temperatures and densities. The Newtonian equations of motion were integrated numerically using the velocity Verlet algorithm with a time step t = 0.001.²¹ Although the mapping of configurations onto local minima are usually defined to emanate from a steepest descent process, this has been found to be numerically very inefficient; instead, we have used a conjugate gradient approach, the first step of which is in fact a steepest descent operation.²² In most cases this more effective procedure identifies the proper inherent structure and is adequate for present purposes. It is worth noting that the Dzugutov pair potential and all its distance derivatives are continuous at the cutoff points a and b, a feature that stabilizes both the Newtonian equation of motion and the conjugate gradient numerical processes; however, this interaction smoothness property has not been present in many past numerical studies that have truncated infinite-range interactions (such as the Lennard-Jones model potential) at a finite-distance cutoff.

All simulations were performed in the microcanonical ensemble, and temperature was adjusted by means of periodic velocity scaling every 10³ time steps. Except for the highest temperature examined (T = 3.0), the system was slowly cooled from reduced temperature T = 2.0 to the desired temperature over the course of 5 × 10⁵ time steps, and then allowed to relax for an additional 2 × 10⁵ steps. During the subsequent production period, velocity rescaling was turned off and energy minimization was performed every 10³ time steps. Twenty time origins, each separated by 10 time units, were used to gather statistics for the calculation of mean-squared displacement curves.

We found that this simulation protocol resulted in the system spontaneously freezing into a somewhat strained and defective crystal at T = 0.8 during the latter portion of the production phase. By using only the data before this freezing event, we were able to extract properties corresponding to the metastable liquid phase at this relatively low temperature. At the slightly lower temperature T = 0.7, however, we found this transition to occur too quickly to extract meaningful liquid-state information. To examine the behavior of the supercooled liquid at even lower temperatures, therefore, we employed a quenching procedure in which an equilibrated high-temperature (T = 1.4) liquid was suddenly forced to a much lower temperature via periodic velocity rescaling. By quenching to low enough temperatures (0.4 T 0.5), crystal nucleation became sluggish owing to the increased viscosity of the melt, permitting adequate measurements to be performed for the liquid state. Visual examination of the final configuration from these trajectories verified that crystallization indeed had not occurred.

Finally, and for comparison with the liquid-phase results, we have investigated the melting behavior of the crystal upon heating at the constant density of = 0.85. By periodic velocity rescaling, we slowly heated a perfect bcc crystal of 250 atoms from T = 0.5 to 1.0 over the course of 10⁶ time steps. Similarly to the liquid-state simulations, inherent-structure minimizations were performed every 10³ time steps, though the mean-squared displacement curves were tabulated for a single time origin.

IV. Results

Figure 1 shows characteristic mean-squared displacement curves for a single trajectory, generated from a system equilibrated at T = 1.0. The discontinuous nature of the inherent-structure curve is immediately obvious; the system remains in its original basin until t = 0.02, after which point it maneuvers rather quickly between nearby basins. The small length scale on which these changes occur suggests that the inherent-structure transitions involve a relatively small number of particles. The corresponding Newtonian calculation remains smooth and exhibits the initial ballistic behavior described by eq II.3 until roughly t = 0.10. It is clear at this temperature that the system leaves its initial basin relatively rapidly, which results in inherent-structure displacements greater than their Newtonian counterparts at subsequent times.

Figure 1 Short-time behavior of a single trajectory for the liquid at = 0.85 and T = 1.0. The Newtonian (solid line) and inherent-structure (line with symbols) trajectories are displayed on a per-particle basis. Each diamond symbol corresponds to a single time step.

Figure 2 illustrates the freezing process which was detected at the slightly lower temperature of T = 0.8. During this trajectory, both the internal energy and pressure decreased abruptly around t = 850, fluctuating around new average values for further molecular dynamics steps. In particular, the pressure drop implies that under isobaric conditions the volume change on melting would be positive. An additional 10³ MD steps did not yield further changes or transitions (not shown). The final inherent-structure configuration is depicted in Figure 3a, which clearly indicates the presence of crystallinity. To collect liquid-state data at lower temperatures and avoid the freezing transition, therefore, a procedure was employed in which high-temperature systems were rapidly quenched and allowed to evolve at T = 0.4 and 0.5 (see preceding section). These temperatures were found to be high enough for the system to eventually escape cagelike dynamics (results not shown), but sufficiently low to suppress crystal nucleation during the trajectory. The final inherent structures from these runs were found to remain amorphous; the T = 0.5 case appears in Figure 3b.

	Figure 2 Progression of the per-particle potential energy and pressure in the T = 0.8 trajectory. At approximately t = 850 the system freezes into a strained and defective crystal.
	Figure 3 Visualizations of the final inherent-structure configurations from (a) the T = 0.8 trajectory, which freezes into a strained and defective crystal, and (b) the T = 0.5 quenched trajectory, which remains in an amorphous state.

We show selected average mean-squared displacement results for both the high-temperature and quench runs in Figure 4. At the relatively high temperature T = 1.6, the vertical position of the inherent-structure displacement curve exceeds that of its Newtonian equivalent by a substantial amount, over 0.1 squared distance units per particle. This again reflects the tendency of the system very quickly to leave its initial basin at increased temperatures, with the minimization procedure on average taking the system to inherent structures that are quite remote. As the temperature is decreased, however, the separation between the two curves gradually decreases and appears to change sign at a temperature that lies between the T = 0.5 and 0.8 cases. At the cooler of these, the vertical ordering of the two curves has reversed. In that case, the system spends a substantial amount of time in each basin before finding a sufficiently low-energy avenue to another, whereas the vibrational displacements around energy minima of the Newtonian trajectory more quickly move the system away from its initial configuration. At all temperatures, the long-time behavior of the two curves yielded statistically identical self-diffusion coefficients, which are shown for the Newtonian version in Figure 5.

	Figure 4 Short-time behavior of the Newtonian (solid line) and inherent-structure (dotted line) trajectories for several temperatures. Each curve is an average over 20 time origins and reported on a per-particle basis. The vertical ordering of the two curves changes between the temperatures T = 0.8 and T = 0.5.
	Figure 5 Arrhenius plot of self-diffusion coefficients from the current work with N = 256 (filled symbols) and that of ref 8 (open symbols) with N = 16000.

A quantitative depiction of the offset between inherent-structure and Newtonian mean-squared displacement results is given in the top of Figure 6. We show calculations averaged from both the entire collection, and from only the first 200 Newtonian/inherent-structure pairs in each trajectory. The two averages reveal a nontrivial amount of statistical fluctuation in the data; however, one can generally extract a decreasing trend in the inherent-structure offset compared to the Newtonian data with decreasing temperature. This relative offset vanishes around T = 0.7 and becomes negative at lower temperatures. In the bottom of this figure, we also show the average squared return distance between instantaneous configurations and their associated inherent structure. The growth of this quantity with temperature is a reflection of the fact that the system spends increasingly more time in basin areas that are remote from inherent structures. Furthermore, its nonlinear behavior is evidently influenced both by vibrational anharmonicity and by changes with temperature of the shapes of occupied basins. Interestingly, it appears this average return distance possesses an inflection point in the same temperature region where the inherent-structure mean-squared displacement offset extrapolates to zero.

Figure 6 Top panel: average difference between the inherent-structure and Newtonian mean-squared displacement curves. The open symbols are averaged over the full trajectory, whereas the closed ones are averaged over only the first 200 time units. Bottom panel: average squared return distance between an instantaneous configuration and its associated inherent structure.

By slowly "heating" a perfect bcc crystal configuration, we were also able to determine the corresponding mean-squared displacement behavior in the solid phase. Presumably, one would expect the inherent-structure displacement to remain zero at low temperatures, where no defects are present in the crystal, and to increase sharply as the temperature rises above the melting point. In simulation, however, the crystal structure is likely to persist as a metastable feature even as the temperature is raised slightly above the true thermodynamic melting line. This is in fact the behavior we observe in our calculations, displayed in Figure 7. The crystal remains stable until the temperature reaches T = 0.95, at which point both the Newtonian and inherent-structure mean-squared displacement curves increase dramatically. (For comparison, the slopes of the liquid-phase curves at that temperature are nearly 5000 times greater than the slope exhibited by the crystal's Newtonian curve prior to melting.) It is interesting to note that melting occurs when <(R - r)²> 0.045, which is equivalent to a Lindemann ratio (of rms displacement to nearest-neighbor separation) equal to 0.18. Although this number may be somewhat inflated as a result of bcc crystal superheating, it is nevertheless within the range reported for other models and substances.^15-17 The average squared return distance in the liquid at T = 1.0 is <(R - r)²> 0.115, which implies a Lindemann ratio approximately equal to 0.29. This is about 1.6 times that of the bcc crystal at the same temperature, thus exhibiting a substantial increase as a result of melting. This behavior is similar to an observation made some years ago for a smoothly truncated version of the Lennard-Jones interaction;¹⁷ however, in that former case the jump factor for the Lindemann ratio was approximately 2.7, implying a significant model dependence. Nevertheless, these results support the applicability of a properly structured "inverse" Lindemann criterion for the liquid phase, in which crystallization spontaneously occurs during liquid cooling when the average return distance of the system to its nearest inherent structure declines to a relevant critical value.

Figure 7 Progression of the Newtonian (solid line) and inherent-structure (dotted line) mean-squared displacement curves as a perfect bcc crystal is heated from T = 0.5 to T = 1.0. When the temperature reaches approximately 0.95, the crystal melts and both displacement curves increase dramatically.

V. Conclusions and Discussion

The long-time limiting behavior of particle mean-square displacements determines self-diffusion constants. Two contrasting routes are available for this determination. The first, and traditional, route monitors the Newtonian dynamical paths for selected particles. The second option relies on mapping the N-particle time-evolving configuration onto inherent structures (potential energy minima), and uses those inherent structures to evaluate mean-square displacements versus time. Although these two alternatives must in principle yield identical results for the self-diffusion constant, they display short and intermediate time differences that arise from details of diffusion processes. These differences hold the promise of revealing insights for deeper understanding of those processes.

The present study involves a series of classical molecular dynamics simulations for the pairwise-additive Dzugutov potential model at reduced density 0.85, and over a wide temperature range that encompasses equilibrium and supercooled liquid states, as well as the crystal. The simulations produced mean-square particle displacement results both for the Newtonian trajectory and for the inherent-structure alternatives. The implied liquid-phase diffusion constants (Figure 5) show a temperature dependence that, within the numerical uncertainty, can be described as an Arrhenius behavior with an activation energy of approximately 3.0 reduced energy units. This activation energy is approximately 5.2 times the depth of the Dzugutov pair potential minimum, eq III.3.

In both protocols, the asymptotic long-time form of the mean-square displacement curves are accurately described as linear functions of time at all temperatures. The best linear fits for the Newtonian and the inherent-structure alternatives exhibit constant offsets C(T) (eq II.4) and C_IS(T) (eq II.8) that, at the conditions investigated here, are positive, but with distinct temperature dependences. The sign of the offset C depends on the relative magnitudes of t₁ = 2Dm/k_BT (the time when 6Dt = <[r(0) - r(t)]²>, assuming parabolic behavior of the mean-squared displacement) and t₂, the momentum relaxation time at which the mean-squared displacement essentially changes from parabolic to linear. If t₂ > t₁, the common situation, C > 0. However, t₂ decreases with increasing density, and hence C can become negative. The manner in which this occurs, and the influence of caging on C, are interesting topics for future investigation.

At low temperatures, in both crystalline and supercooled liquid states, C > C_IS, due primarily to inclusion of intrabasin displacements in C but their exclusion from C_IS. At high temperature the order is reversed, C_IS > C. In this latter regime, typical liquid configurations are widely removed from their mapped inherent structures on account of strong local density and structural fluctuations driven by intense thermal excitation; by contrast the Newtonian trajectories are immune to this effect. Our results suggest that the crossover temperature at which C = C_IS is approximately T = 0.7. This crossover is analogous to, but not necessarily to be directly identified with, the mode-coupling-theory structural arrest temperature T_x, below which only activated hops between neighboring basins are available for relaxation and diffusion.²³ In fact, Dzugutov and co-workers⁸ have calculated T_x to be about 0.39, and their results seem to suggest that the C(T), C_IS(T) crossover better coincides with the onset of stretched-exponential relaxation and the formation of large icosahedral clusters.

A natural byproduct of the numerical construction of inherent-structure mean-square displacements versus time is evaluation of the intrabasin root-mean-square magnitude of particle return distance under the mapping to potential minima. In the crystal phase, the dimensionless ratio of this length to the nearest-neighbor lattice spacing plays a fundamental role in the Lindemann melting criterion.^15,16 For the Dzugutov potential model at density 0.85 we have verified that the Lindemann ratio for the stable body-centered cubic crystal rises from zero at absolute zero temperature to approximately 0.18 at the melting point, in agreement with conventional expectation. Furthermore, the Lindemann ratio in the liquid at the freezing point is substantially larger. This latter observation is similar to one made some time ago for a modified Lennard-Jones potential model,¹⁷ and thus strengthens the case for advocating a suitably crafted "inverse-Lindemann" freezing criterion. A significant jump in the Lindemann ratio upon passing from crystal to liquid at the melting point accents the qualitative difference in character of the potential energy basins that respectively underlie the two phases.

Finally, we remark briefly on the comparison between the present results for the Dzugutov model at reduced density 0.85 to the Keyes and Chowdhary simulation results (N = 32) for the Lennard-Jones model at reduced density 1.0.⁹ First, those authors have also observed distinct offsets for Newtonian and for inherent-structure mean-square displacement curves as a function of time, presenting qualitatively the same kinds of temperature variations, with a crossover temperature of approximately 0.8 (Figure 10 of ref 9). Keeping in mind the relative depths of the Lennard-Jones and Dzugutov pair potentials, however, this temperature is effectively lower than that of the present Dzugutov study. Second, the self-diffusion rates for that Lennard-Jones study (Figure 7 of ref 9) appear to have a range of activation energies that are a considerably smaller multiple (roughly 2.0) of the pair potential depth than we have observed for the Dzugutov interaction model. This may seem somewhat surprising, given the higher packing density of the Lennard-Jones case (even after accounting for the slightly different distances at which the two potentials have first zeros, 1.0000 for Lennard-Jones and 1.0245 for Dzugutov). In addition, the Dzugutov model may have a greater propensity for heterogeneous dynamics.⁸ Evidently, the details of short-range order that are distinctly different in the respective liquids, stemming from the contrasting interactions, markedly influence the diffusion processes in the two cases. In short, these two models present a clear and instructive contrast.

Acknowledgment

P.G.D. gratefully acknowledges financial support by the U.S. Department of Energy, Division of Chemical Sciences, Geosciences, and Biosciences, Office of Basic Energy Sciences, Grant No. DE-FG02-87ER13714. M.S.S gratefully acknowledges the support of the Fannie and John Hertz Foundation.