An \emph{ab initio} method for locating characteristic potential energy minima of liquids
Abstract
It is possible in principle to probe the many–atom potential surface using density functional theory (DFT). This will allow us to apply DFT to the Hamiltonian formulation of atomic motion in monatomic liquids [Phys. Rev. E 56, 4179 (1997)]. For a monatomic system, analysis of the potential surface is facilitated by the random and symmetric classification of potential energy valleys. Since the random valleys are numerically dominant and uniform in their macroscopic potential properties, only a few quenches are necessary to establish these properties. Here we describe an efficient technique for doing this. Quenches are done from easily generated “stochastic” configurations, in which the nuclei are distributed uniformly within a constraint limiting the closeness of approach. For metallic Na with atomic pair potential interactions, it is shown that quenches from stochastic configurations and quenches from equilibrium liquid Molecular Dynamics (MD) configurations produce statistically identical distributions of the structural potential energy. Again for metallic Na, it is shown that DFT quenches from stochastic configurations provide the parameters which calibrate the Hamiltonian. A statistical mechanical analysis shows how the underlying potential properties can be extracted from the distributions found in quenches from stochastic configurations.
pacs:
05.70.Ce, 61.20.Gy, 71.15.NcI Introduction
The potential energy surface underlying the motion of a monatomic liquid is composed of intersecting valleys in the many–atom configuration space. In Vibration–Transit (VT) theory, the atomic motion is described by vibrations within a single valley, interspersed by transits, which carry the system across intervalley intersections Wallace (1997); Chisolm and Wallace (2001). The pure vibrational motion is described by a tractable Hamiltonian which accounts for the dominant part of the liquid’s thermodynamic properties. The remaining transit contribution is complicated but small. The vibrational Hamiltonian is calibrated by potential parameters evaluated at local minima (called structures) of the potential surface. These have been previously evaluated for a model of liquid Na based on an interatomic pair potential De LorenziVenneri and Wallace (2007). Our goal now is to introduce first principles electronic structure calculations within density functional theory (DFT) to calculate the structural and vibrational parameters of the VT Hamiltonian.
In recent years, DFT has been used successfully in liquid studies to support the development of ab initio Molecular Dynamics (MD). The original Car–Parrinello method Car and Parrinello (1985) was successfully applied to the melting of Si Sugino and Car (1995). Calculations of the MD trajectory on the adiabatic (electronic ground state) potential surface have been performed for a broad spectrum of elemental liquids Kresse and Hafner (1993a, b); Kresse (2002). This work has provided accurate and physically revealing results for Al de Wijs et al. (1998), Fe Alfè et al. (2000), and Ge Chai et al. (2003). The techniques we employ to calculate the supercell ground state energy and Hellmann–Feynman forces are quite similar to those of Kresse et al. Kresse and Furthmüller (1996). On the other hand, rather than proceed to MD calculations, our objective is to calibrate the (dominant) vibrational part of the liquid dynamics Hamiltonian from properties of local potential minima. We believe that the Hamiltonian formulation will usefully complement ab initio MD in the study of dynamical properties of liquids.
The procedure of quenching a system to its potential energy minima was introduced by Stillinger and Weber Stillinger and Weber (1982); Weber and Stillinger (1984); Stillinger and Weber (1984), and has become a valuable technique for studying systems with interatomic potentials. The traditional method is to quench to many structures from an equilibrium MD trajectory, and then use statistical mechanics to extract the underlying statistical properties of the potential surface Heuer and Büchner (2000); Sciortino et al. (1999). In VT theory, however, we need very few structures—in principle only one for a given liquid at a given volume. Hence the traditional procedure for finding structures is inefficient, as a large number of MD iterations is required to bring the system initially to equilibrium as well as to avoid unwanted correlations between quenches. The problem is especially severe for DFT, where each iteration requires a converged total energy calculation, which is computationally very costly.
We propose a simpler and more efficient method to probe the underlying potential landscape. Rather than quench from equilibrium MD configurations, we quench from configurations that are independent of interatomic interactions, and very fast to generate. Our purpose here is to demonstrate two properties of this quench method: that it is capable of producing the entire distribution of potential energy minima, and when used with DFT it can achieve our goal of ab initio calibration of the Hamiltonian.
To perform the calibration, we will interpret structural data via the original hypothesis of VT theory Wallace (1997): The potential energy valleys are classified as random and symmetric. The random valleys numerically dominate the liquid statistical mechanics as , and they all have the same potential energy properties as ; hence any one such valley may be used to calibrate the Hamiltonian. The symmetric valleys have a broad range of potential energy properties, but make an insignificant contribution to the liquid statistical mechanics as . This hypothesis has been verified for the pair potential model of liquid Na at De LorenziVenneri and Wallace (2007), and work in progress extends this verification to Holmström et al. (2009).
The quench calculations are carried out for metallic Na at the density of the liquid at melt, using our model interatomic potential, and also using DFT. The Na potential was derived in pseudopotential perturbation theory, with an added Born–Mayer repulsion, and was calibrated from experimental crystal data at zero temperature and pressure Wallace (1968). The potential has since been shown to give excellent results for a broad range of experimental properties of crystal and liquid Na (for a partial summary, see Wallace (2002)). While there is no doubt DFT will provide accurate total energy results, we shall still need to verify that DFT quenches arrive at random structures, not symmetric ones. This verification will be accomplished with the aid of the Na interatomic potential results.
Our application of DFT to liquid dynamics theory is being pursued for a number of nearly–free–electron metals and transition metals. A preliminary report has been presented on results for Na and Cu Bock et al. (2008). We have not attempted to study elemental liquids whose equilibrium configurations are influenced by molecular bonding. Examples are As, Se, and Te, which have strong and weak bonds, and Ge, whose liquid structure shows a contribution from covalent crystal bonds. This anisotropic bonding will complicate the random structures underlying the motion in such liquids. This complication remains beyond the scope of the present work. Moreover, since we have not yet presented structural data for other liquid metals, the present conclusions are strictly valid only for liquid Na. The results are expected to apply to many elemental liquids, perhaps all of them, but this extension is not demonstrated here.
We consider a system of atoms in a cubic box with periodic boundary conditions. We construct configurations in which the nuclei are distributed uniformly over the box, within a constraint limiting the closeness of approach of any pair. These are called constrained stochastic, or simply stochastic, configurations. The procedure is described in Sec. II, and the spatial uniformity is verified by means of pair distribution functions. In Sec. III, our twofold purpose is addressed by two separate quench studies. In Sec. III.1, the Na pair potential is used to quench both equilibrium MD configurations and stochastic configurations. Comparison of the results will validate the use of stochastic configurations. In Sec. III.2, DFT is used to quench Na from stochastic configurations. Comparison of the results with pair potential results will confirm that the DFT structures are random and therefore calibrate the Hamiltonian. In Sec. IV, relations are derived between quenched distributions and the densities of states in the underlying potential energy surface. This analysis provides the statistical mechanical framework for interpreting results of the present quench technique. Our conclusions are summarized in Sec. V.
Ii Generating stochastic configurations
Only minimal information is required to generate stochastic configurations: the number of atoms and the system volume , plus a distance of closest approach which is described below. Nothing of boundary conditions or the system potential has to be specified; however, after the stochastic configuration is constructed, for all further calculations periodic boundary conditions are used.
For a cubic cell with volume , we construct a configuration by choosing the particle coordinates at random over the cell. Randomness is important, as we shall use it in Sec. IV to determine the statistical weight factors for stochastic configurations. Next, a configuration is discarded if any two atoms are closer than a distance . This is done for practical reasons: The self consistent field (SCF) calculation of DFT will not converge if atoms are too close to each other, and the pair potential at very small radii could lead to numerical instabilities in the conjugate gradient method due to the repulsive core. In practice, the excluded space can be very small. For Na we choose Å, so the relative excluded space , where is the volume per atom, is only . Hence the stochastic configurations are expected to be spatially uniform to a very good approximation.
To test the uniformity of stochastic configurations, we examine their pair distribution functions . The conditional probability density is constructed as follows: Pick a system atom as central atom and denote its position by . Make a set of bins labeled , in the form of concentric shells. Bin has inner radius , outer radius , and volume . The pair distribution function has histogram , where is the number of atoms in bin . Given the small size of described in the previous paragraph, we expect to be nonzero even at relatively small radii. It is therefore important to normalize the bin count of bin with the correct volume of the bin, instead of using the approximate volume as is often done. The bin contents are then averaged over the choice of each system atom as the central atom. While the bin radii are arbitrary, we usually take either constant or constant.
For Na at , we constructed 1000 stochastic configurations and the histogram for each. With constant, the mean and the standard deviation of the histogram are shown in Fig. 1. The blank space at small is the empty sphere of radius . The scatter at small reflects the decreasing as decreases, and the corresponding decrease in . With constant, the mean and standard deviation of the histogram is shown in Fig. 2. There the distance between points increases as decreases, but the standard deviation remains nearly constant because remains nearly constant. The figures show the uniformity of stochastic configurations for , with being very small compared to the mean nearest neighbor distance.
Iii Properties of the quenched structures
iii.1 Validation of quenches from stochastic configurations
Figs. 4 and 5 compare two distributions of , each from 1000 pair potential quenches at . Fig. 4 is obtained by steepest–descent quenches from equilibrium MD at 800 K. This figure is an extension of the work reported by De LorenziVenneri and Wallace (2007). Fig. 5 is obtained by conjugate gradient quenches from stochastic configurations. We have verified the equivalence of the two quench techniques for our system; for a related verification, see Chakravarty et al. (2005); Chakraborty and Chakravarty (2006).
In each histogram we see the distinct random and symmetric distributions, consistent with the VT hypothesis. The random distribution is taken to be the dominant peak, out to where the histogram vanishes on either side. The symmetric structures are interpreted as the isolated parts of the histogram that are located outside the main peak on the low–energy side.
The random distribution is numerically dominant and very narrow. The mean and standard deviation of each random distribution is given by
(1) 
Note that the dominant volume–dependent part of the pair potential is omitted here (see Wallace and Clements (1999), Eq. (1.1) and Fig. 1). For this reason, most of the binding energy is missing from the pair potential energies in Eq. (1). The mean value is the most accurate approximation to the thermodynamic limit value that we can make. The standard deviation is the error expected from quenching only once and using that result. If is increased, the mean value is expected to change slightly in converging to its thermodynamic limit, while the standard deviation goes to zero as . For a physical measure of the difference in mean values, we note that the main contribution to the liquid thermal energy is the classical vibrational contribution , which is 95.90 meV/atom at . The difference in means is 0.08% of this. Experimental error in the thermal energy of elemental liquids at is typically (0.1  0.5)%. Hence the two random distributions in Figs. 4 and 5 are identical to better than experimental error.
Performing so many quenches has allowed us for the first time to see a clear and meaningful distribution of symmetric structures (compare for example De LorenziVenneri and Wallace (2007); Clements and Wallace (1999); Wallace and Clements (1999)). Quenching from equilibrium MD yields 18 symmetric structures out of 1000 quenches, Fig. 4. Quenching from stochastic configurations yields 23 symmetric structures in 1000 quenches, Fig. 5. Very approximately, the symmetric distribution is constant and ranges from the bcc crystal, meV/atom Wallace and Clements (1999), to the lower end of the random distribution. This broad distribution with few structures is consistent with the VT hypothesis. If is increased the symmetric distribution width is expected to remain the same, while the relative number of symmetrics is expected to become negligible. In all these properties, the symmetric distributions in Figs. 4 and 5 are the same to statistical accuracy.
iii.2 Calibration of the VT Hamiltonian
In order to calculate ab initio the thermodynamic properties of liquid Na for comparison with experiment, we quenched a stochastic configuration at Å with DFT Bock et al. (2008). The normal mode frequency spectrum was also calculated at this volume. For almost all monatomic liquids, the vibrational motion is nearly classical at . This means that the essential information required from is the logarithmic moment of the frequencies, which provides the characteristic temperature . Hence for calculation of thermodynamic properties, the VT Hamiltonian is calibrated from for the energy, and for the entropy. Additional data which automatically accompanies the calculation of will calibrate the Hamiltonian for nonequilibrium statistical mechanics.
The DFT calculation is done with the VASP code VAS (), using the projector augmented wave (PAW) method in the generalized gradient approximation (GGA) Blöchl (1994); Kresse and Joubert (1999). The planewave energy cutoff is 101.7 eV, the maximum core radius is 2.5 Å, and the total energy convergence criterion is eV. We use a –centered Monkhorst–Pack grid Monkhorst and Pack (1976) with 14 –points in the irreducible Brillouin zone. The total energy convergence for these parameters was carefully verified. Note that it is the large size of the real–space supercell which allows us to use few –points in comparison to the large number (several thousands) needed for crystal metal calculations with small unit cells. The quench is done by nonlinear conjugate gradient method. The system is considered quenched when the energy difference between subsequent configurations is 10 eV or less. The DFT structure is at , a number large enough to get potential energy parameters accurate to a few percent, but small enough that convergence properties of the calculations can be studied.
Because of the strong dominance of random valleys in the potential energy surface, the DFT structure is expected to be random. To eliminate different zeros of energy, we evaluate the energy difference
(2) 
where the superscripts and represent respectively the random structure and the bcc crystal. The comparison is
(3) 
The value for the pair potential is from the second mean in Eq. (1), which is also calculated from stochastic configurations. The difference of meV/atom is small compared to theoretical errors in , and also compared to experimental error in the energy of liquid Na at melt.
An independent confirmation is furnished by . The Na pair potential value is from De LorenziVenneri and Wallace (2007):
(4) 
The relative difference of in will make a corresponding difference of in the theoretical entropy of liquid Na at melt. The difference is well within theoretical error in , and is close to the experimental error in the entropy.
The structural pair distribution is not a parameter of the VT Hamiltonian, but has a role in density fluctuation phenomena, and it is therefore interesting to compare the DFT and pair potential results. The comparison is shown in Fig. 6, where the agreement is excellent. Notice the DFT curve () has a small deficiency at the tip of the first peak, compared to the pair potential curve (). This deficiency is a small– effect, and is observed also with the pair potential at , but not at (Clements and Wallace (1999), Fig. 2).
Iv Extracting densities of states from quench results
iv.1 Quenches from equilibrium configurations
In classical statistical mechanics, the partition function for a single potential valley harmonically extended to infinity is . The factor expresses the vibrational motion. The transit contribution, which accounts for the valley–valley intersections, will be neglected here. The total liquid partition function is
(5) 
where is the joint density of states for the collection of valleys. The normalization of is , the total number of valleys. The equilibrium statistical weight of a single valley is
(6) 
In equilibrium at , the probability of finding the system in at , and in at , is , where
(7) 
Upon quenching from an equilibrium trajectory at , the structures sampled will exhibit a distribution proportional to foo (). In view of Eqs. (5) and (6) it follows that is insensitive to the normalization of . Therefore measurements of can not be used to count the valleys.
The probability of finding the system in at is , where
(8) 
Upon quenching from an equilibrium trajectory at , the structures sampled will exhibit a distribution of proportional to . The distribution in Fig. 4 is proportional to . However, the density of states in is , given by
(9) 
and differs from by the statistical weight . In our studies of liquid dynamics, the purpose of quenching is to determine the densities of states and , because these are parameters of the Hamiltonian. These densities of states must be solved for from the observed distributions and from the above equations. Even though the symmetric structures are supposed to be unimportant for the liquid as , those structures will always be statistically important at , and will therefore be included in our analysis.
Let us introduce subscripts and for random and symmetric, respectively, and write
(10) 
and correspondingly . From Figs. 4 and 5, at has very small width, comparable to experimental error in the internal energy of the liquid at melt, and this width is expected to decrease further as increases Holmström et al. (2009). These results suggest a model for . Let us define the liquid as the mean for random structures when , with a similar definition for . The model is
(11) 
where , and is the number of random valleys. From this it follows that , and the random contributions to Eqs. (7) and (8) become
(12) 
Hence, to finite– errors, the random valley Hamiltonian parameters and are determined directly by quenches from equilibrium MD at . And, because of the form of Eq. (11) for , these observations will remain true when the statistical mechanics theory is improved to include transit effects.
The symmetric density apparently has –independent width with ranging from to . Symmetric structures with exist, but they are rare for monatomic systems. The dependence of is not trivial. One expects that additional (symmetry) parameters are important for symmetric structures. Nevertheless, and are well defined, and can be extracted from quench data with the aid of Eqs. (7) and (8).
iv.2 Quenches from stochastic configurations
On an equilibrium trajectory at , the probability the system is found in a given potential energy valley is for that valley. The statistical weight is quite different for stochastic configurations. These configurations are uniformly distributed over configuration space, except for the small excluded Cartesian–space volume at each nucleus. Neglecting this constraint, the probability the system is found in a given potential valley is the valley volume divided by the entire –dimensional volume. Let us denote the corresponding statistical weight factors as and for random and symmetric valleys respectively.
Because of their uniformity, the random valleys all have the same volume in the thermodynamic limit. To arrive at a complete solution, it is necessary to include the symmetric valleys. Let us assume that they also have a uniform configuration space volume. The number of random (symmetric) valleys is denoted (), and the single–valley volume is (). The statistical single–valley weights are
(13)  
(14) 
The probability distributions are
(15)  
(16) 
Hence the random and symmetric densities of states are each proportional to the distribution found in quenches from stochastic configurations. Applying the model of Eq. (11) for yields
(17) 
The conclusion here is the same as with equilibrium configurations, Eq. (12), that the parameters and are determined directly by quenches from stochastic configurations, up to finite– errors. The reason, of course, is the form of the model for , Eq. (11). For symmetric structures, the above equations reveal two significant points:

Quenches from stochastic configurations can determine the magnitude of relative to , but only when is known.

The relation between as determined by the two quench methods is unknown until is evaluated.
These points are relevant to the distributions shown in Figs. 4 and 5.
V Conclusions
In this article, we have addressed two main research goals: (1) The use of stochastic configurations to probe the distribution of potential energy minima, and (2) the calibration of the VT Hamiltonian with DFT.
v.1 Stochastic Configurations
Quenches from stochastic configurations produce statistically indistinguishable distributions of the potential energy, , compared to quenches from equilibrium liquid MD trajectories. We have demonstrated that quenching from stochastic configurations can be used to find the entire distribution of Na potential energy minima, i.e. they reliably reproduce the correct distribution of random and symmetric structures. This is illustrated by a comparison of distributions for quenches from equilibrium MD and from stochastic configurations (Figs. 4 and 5, Sec. III.1). Quenches from stochastic configurations yield an accurate random distribution for Na in agreement with the random distribution from quenches from equilibrium MD, Eq. (1).
Compared to generating and selecting configurations from equilibrium MD trajectories, our stochastic configuration method is significantly faster and more economical (Sec. II). Stochastic configurations do not require interatomic potentials or costly equilibration and long simulation times. Simply generate random Cartesian coordinates for each atom under a minimal excluded–volume constraint to eliminate particle overlap. Hence, this procedure requires very little computational effort, an economy that accommodates ab initio quench methods even for large systems.
v.2 Calibration of the VT Hamiltonian
Calibration of the VT Hamiltonian is based on the presumed dominance and uniformity of random valleys as . This view is given mathematical expression in the model for , Eq. (11). It follows that the thermodynamic limit parameters and are determined directly from data for either MD quenches or stochastic quenches, up to finite– errors.
We have demonstrated that the DFT structure in Sec. III.2 is random by comparing the mean potential energy with the pair potential results in Eq. (3), the phonon moment in Eq. (4), and the pair distribution function in Fig. 6. We conclude that being random, the DFT structure can provide ab initio calibration of the VT Hamiltonian.
As verified by their pair distribution function, stochastic configurations have nuclei distributed nearly uniformly over the system volume (Figs. 1 and 2, Sec. II). Therefore among stochastic configurations, the statistical weight for a many–atom potential energy valley is (nearly) proportional to the valley volume. This is in contrast to quenches from equilibrium MD, which require extensive modeling to extract the Boltzmann factor from the weight Heuer and Büchner (2000). Given the statistical weight, characteristics of the underlying potential surface can be extracted from data acquired by stochastic quenches (Sec. IV).
We do not suggest that DFT quenches from stochastic configurations will invariably arrive at random structures, just as quenches from an equilibrium MD trajectory may result in a symmetric structure. Indeed, some symmetric structures have appeared in our DFT quenches. Precisely what is required to eliminate symmetric structures from any collection of quench data is an ongoing research question.
Vi Acknowledgments
EH and RL would like to thank Fondecyt projects 11070115 and 11080259. They are also thankful for the support of “Proyecto Anillo” ACT24. This work was funded by the US Department of Energy under contract number DEAC5206NA25396. The open and friendly scientific atmosphere of the Ten Bar International Science Café is also acknowledged.
Footnotes
 preprint: LAUR 0901331, LAUR 0903784
References
 D. C. Wallace, Phys. Rev. E 56, 4179 (1997).
 E. D. Chisolm and D. C. Wallace, J. Phys.: Condens. Matter 13, R739 (2001).
 G. De LorenziVenneri and D. C. Wallace, Phys. Rev. E 76, 041203 (2007).
 R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985).
 O. Sugino and R. Car, Phys. Rev. Lett. 74, 1823 (1995).
 G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993a).
 G. Kresse and J. Hafner, Phys. Rev. B 48, 13115 (1993b).
 G. Kresse, J. NonCryst. Solids 312–314, 52 (2002).
 G. A. de Wijs, G. Kresse, and M. J. Gillan, Phys. Rev. B 57, 8223 (1998).
 D. Alfè, G. Kresse, and M. J. Gillan, Phys. Rev. B 61, 132 (2000).
 J.D. Chai, D. Stroud, J. Hafner, and G. Kresse, Phys. Rev. B 67, 104205 (2003).
 G. Kresse and J. Furthmüller, Computational Materials Science 6, 15 (1996).
 F. H. Stillinger and T. A. Weber, Phys. Rev. A 25, 978 (1982).
 T. A. Weber and F. H. Stillinger, J. Chem. Phys. 80, 2742 (1984).
 F. H. Stillinger and T. A. Weber, Science 225, 983 (1984).
 A. Heuer and S. Büchner, J. Phys.: Condens. Matter 12, 6535 (2000).
 F. Sciortino, W. Kob, and P. Tartaglia, Phys. Rev. Lett. 83, 3214 (1999).
 E. Holmström, N. Bock, T. B. Peery, R. Lizárraga, G. De LorenziVenneri, E. D. Chisolm, and D. C. Wallace (2009), work in progress.
 D. C. Wallace, Phys. Rev. 176, 832 (1968).
 D. C. Wallace, Statistical Physics of Crystals and Liquids (World Scientific, New Jersey, 2002).
 N. Bock, T. Peery, E. D. Chisolm, G. De LorenziVenneri, D. C. Wallace, E. Holmström, and R. Lizárraga, Bull. Am. Phys. Soc. 53 (2), J9:00004 (2008).
 C. Chakravarty, P. G. Debenedetti, and F. H. Stillinger, J. Chem. Phys. 123, 206101 (2005).
 S. N. Chakraborty and C. Chakravarty, J. Chem. Phys. 124, 014507 (2006).
 D. C. Wallace and B. E. Clements, Phys. Rev. E 59, 2942 (1999).
 B. E. Clements and D. C. Wallace, Phys. Rev. E 59, 2955 (1999).
 Vienna Abinitio Simulation Package (VASP), URL http://cms.mpi.univie.ac.at/vasp/.
 P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
 G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
 H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
 An equilibrium MD system moves at constant total energy and total momentum. These constraints affect fluctuations in order 1, but affect mean values only in order Lebowitz et al. (1967). We exercise caution when using canonical weighting to analyze MD results.
 J. L. Lebowitz, J. K. Percus, and L. Verlet, Phys. Rev. 153, 250 (1967).