Fast Non-Adiabatic Dynamics of Many-Body Quantum Systems
Brett Larder, Dirk Gericke, Scott Richardson, Paul Mabey, Thomas White, Gianluca Gregori
FFast Non-Adiabatic Dynamics of Many-Body Quantum Systems
B. Larder, D.O. Gericke, S. Richardson,
1, 3
P. Mabey, T. G. White, and G. Gregori ∗ Department of Physics, University of Oxford, Parks Road, Oxford OX1 3PU, UK Centre for Fusion, Space and Astrophysics, Department of Physics,University of Warwick, Coventry CV4 7AL, UK AWE, Aldermaston, Reading, Berkshire RG7 4PR, UK LULI – CNRS, Ecole Polytechnique, CEA, Universit´e Paris-Saclay, F-91128 Palaiseau Cedex, France Department of Physics, University of Nevada, Reno, Nevada 89557, USA
Modeling many-body quantum systems withstrong interactions is one of the core challengesof modern physics. A range of methods has beendeveloped to approach this task, each with itsown idiosyncrasies, approximations, and realm ofapplicability[1]. Perhaps the most successful andubiquitous of these approaches is density func-tional theory (DFT)[2]. Its Kohn-Sham formula-tion [3] has been the basis for many fundamen-tal physical insights, and it has been successfullyapplied to fields as diverse as quantum chem-istry, condensed matter and dense plasmas [4–10]. Despite the progress made by DFT and re-lated schemes, however, there remain many prob-lems that are intractable for existing methods.In particular, many approaches face a huge com-putational barrier when modeling large numbersof coupled electrons and ions at finite tempera-ture [10–12]. Here, we address this shortfall witha new approach to modeling many-body quan-tum systems. Based on the Bohmian trajecto-ries formalism [13, 14], our new method treatsthe full particle dynamics with a considerable in-crease in computational speed. As a result, weare able to perform large-scale simulations of cou-pled electron-ion systems without employing theadiabatic Born-Oppenheimer approximation.
Let us consider a many-particle electron-ion systemat finite temperature. In calculating the dynamics ofboth the electrons and ions, we must account for the factthat the ions evolve multiple orders of magnitude moreslowly than the electrons, as a result of their much highermasses. If we are interested in the long-time ionic dy-namics (for example, the ion mode structure), we face achoice of how to deal with this timescale issue. We can ei-ther model the system on the time scale of the electrons– non-adiabatically – and incur a significant computa-tional cost (a cost that is prohibitive in most simulationschemes); or we can model the system on the timescaleof the ions – adiabatically – by treating the electrons as astatic, instantaneously adjusting background. The latterapproach is far cheaper computationally, but does notallow for a viable description of the interplay of ion andelectron dynamics. ∗ Corresponding author: [email protected]
The method we propose here enables us to employthe former (non-adiabatic) approach, retaining the dy-namic coupling between electrons and ions, by signifi-cantly reducing the simulation’s computational demands.We achieve this by treating the system dynamics withlinearized Bohmian trajectories. Having numerical prop-erties similar to those of molecular dynamics for classicalparticles, our approach permits calculations previouslyout of reach: systems containing thousands of particlescan be modeled for long (ionic) time periods, so that dy-namic ion modes can be calculated without discountingelectron dynamics.To begin constructing our method, we consider Bohm’sformulation of quantum mechanics. We can imaginea classical N -body system as a “ N -trajectory” movingthrough 3 N -dimensional configuration space. Bohm’stheory treats an N -body quantum system as an ensembleof these classical N -trajectories, interacting through anadditional N -body potential V B . This Bohm potential isa functional of the density of N -trajectories in configu-ration space, Φ, and a function of the spatial position x ;that is, V B = V B ( x | Φ). Provided matching initial con-ditions, the Bohm ensemble of N -trajectories reproducesthe dynamics of the probability density – as given by theSchr¨odinger equation – exactly [13] (see also Supplemen-tary Information for more details).In its exact form, Bohm’s formulation is as intractableas the N -body Schr¨odinger equation, requiring simula-tion of a huge number of N -body interacting classicalsystems. However, we can construct a fast computa-tional method solving Bohm’s theory by introducing athermally averaged, linearized Bohm potential. The ex-act (but inaccessible) calculation for a pure quantumstate with many particles – based on the theory above –would require us to propagate an array of N -trajectoriesthrough time, at each step recalculating their density andthereby V B (see Fig. 1a). Reliable calculations of den-sity in 3 N -dimensional space would require a prohibitivenumber of N -trajectories, however, making this unfeasi-ble. Here, we propose an alternative: as opposed to ap-plying this theory to a pure quantum state, we considerpropagating an array of thermally coupled N -trajectoriesin a similar manner, as a model of a system at finite tem-perature. Our core assumption is that the time evolutionof a finite-temperature quantum system can be approx-imated by a similar procedure to that used for the purestate; we consider an array of N -trajectories, each cou- a r X i v : . [ c ond - m a t . s t a t - m ec h ] N ov pled to a heat bath setting its temperature, evolving un-der a potential that is a functional of N -trajectory den-sity in configuration space. This procedure can be seenas a trajectory-based analogue of the linearization of theBohm potential over states in quantum hydrodynamics[16].For simplicity, we focus initially on systems in ther-mal equilibrium. We allow our thermal N -trajectories– together modeling the probability density of our finite-temperature system – to evolve in time under a linearizedmean Bohm potential of the underlying pure states. Inthis linear approximation, we replace the mean Bohm po-tential experienced, expressible as a sum over functionalsof individual states, with a functional of a sum over in-dividual states.In addition to moving focus from a pure state to amore practically important finite-temperature state, thekey feature of our approximation is that it dramaticallyreduces the computational expense required to simulatethe system (as compared to the exact pure state case).Crucially, now that we are considering finite-temperature N -trajectories, we need only track a single N -trajectorythrough time, rather than an impractical number ofthem. This follows from two observations:
1. Properties of a classical system with correlation-dependent potentials can be determined self-consistently.
For an arbitrary system of N well-localized particles,the N -particle correlation function can be written as g ( x ) = P ( x ) /P ( x ), where P denotes the joint positionalprobability distribution of the particles, and P is thedistribution for a non-interacting classical system withequal particle densities. Here, the variable x is the setof particle positions, x = { x , x , ..., x N } . We may con-struct inter-particle potentials that are functionals of g : V ( x ) = V ( x ) + V g ( x | g ), where V denotes pair interac-tions and external forces, and V g is a contribution thatvaries with g . The equilibrium properties of a systemin this potential can be found self-consistently. Start-ing from an initial guess for V ( x ), we can calculate g ( x ),through a Monte Carlo simulation or a similar method[15]. This value of the N -particle correlation function g then gives rise to a new approximation for the potential.Iterating this procedure allows for both g and V to befound.
2. Our linearized Bohmian system is equivalent toa classical system with correlation-dependent potentials.
Consider a number of coupled thermal N -trajectoriesrepresenting our linearized Bohmian system. Assumingthat the system is in a temperature regime in which theparticle motion is ergodic, we find that each N -trajectoryhas the same time-integrated correlations; that is, eachhas the same g . In the limit of infinite N -trajectoriesin our ensemble, it follows that the configuration spacedensity of N -trajectories Φ is exactly proportional to thiscommon g . As a result, each N -trajectory moves in acommon static potential (see Fig. 1b), and, as this staticpotential is a functional of configuration space density,Φ, it is equivalently just a functional of g . FIG. 1. Schematic of the applied linearization approxima-tion. a) The time evolution of an N -trajectory in an exactBohmian representation of a pure quantum state (upper), andthe Bohm potential, V B , that it experiences (lower). V B is afunctional of the density of N -trajectories in configurationspace, Φ. At each time step, all of the N -trajectories in theensemble must be updated, Φ calculated, and the updatedBohm potential determined. b) The time evolution of an N -trajectory in our linearized Bohmian representation of a ther-mal state (upper), and the Bohm potential that it experiences(lower). We need only track a single N -trajectory: its cou-pling to a heat bath ensures its ergodicity, so that Φ becomesequal, in equilibrium, to its time-averaged density in config-uration space. As a result, the N -trajectory evolves with atime-independent Bohm potential, generated self-consistentlyby its own time-integrated density. c) The block panel sum-marizes our scheme to determine the Bohmian dynamics andgives the sources for the different potentials needed as input. When combining the results above, our linearizationapproximation becomes a simple mapping: V B ( x | Φ) V g ( x | g ) = − λ ~ √ g N X i =1 ∇ i m i √ g , (1)where m i is the mass of particle i and λ is a lineariza-tion factor that still needs to be determined. As g iscommon to all N -trajectories, our approximation schemeallows us to consider just a single N -body classical sys-tem (Fig. 1b). The required simulation is thus amenableto (computationally cheap) classical Molecular Dynam-ics [15]. The classical particle trajectories simulated thenapproximate the statistics of the full quantum system.Before the scheme above can be implemented, we mustovercome a final fundamental hurdle: the full correlationfunction g appearing in Eq. (1) is too complicated tobe modeled directly (similar to the full N -body wave-function). Therefore, we seek an approximate closure forthis object in terms of lower-order correlations, whichcan be calculated accurately. For this goal, we employthe pair product approximation, whereby the N -bodycorrelation is replaced by a product of pair correlations.Further, we generalize the dependence on λ to a set of λ ij to accommodate different particle species in the pairinteractions.We need an additional correction to our potentials tofulfil the spin statistics theorem, as the Schr¨odinger equa-tion – and thus Bohmian mechanics – does not incorpo-rate particle spin directly. In particular, this correctionwill generate a Fermi distribution for the electrons inthermal equilibrium. Similar to successful approaches ap-plied in quantum hydrodynamics and classical map meth-ods [16–18], we introduce an additional Pauli potentialterm. This term is constructed such that exchange ef-fects are reproduced exactly for a reference electron gassystem. We also employ pseudo-potentials, commonlyapplied in modern DFT calculations, to represent coreelectrons bound in deep shells of the ions by an effectiveion potential seen by the valence electrons.Finally, it remains to set the linearization parameters λ ij to fully determine the system’s Hamiltonian. In thiswork, we take λ ij = 1 for the ion-ion and ion-electronterms. To determine the electron-electron parameter λ ee ,we match static ion correlations – that is, pair distri-bution functions – obtained by DFT calculations. Thismatching can be carried out rigorously with a general-ized form of inverse Monte Carlo (see Supplement Infor-mation). In this way, we determine λ ee with only static information of the system. Subsequent dynamic simula-tions can then be carried out without any free parame-ters.We can now implement the Bohmian dynamics methodwith a Molecular Dynamics simulation inclusive of thepotentials discussed above. Within the microcanonicalensemble, we could simply integrate the equation of mo-tion for the particles. However, we want to simulate sys-tems with a given temperature which requires the cou-pling of the system to a heat bath. Here, we employa modified version of the Nos´e-Hoover thermostat. Itsstandard form is popular in classical Molecular Dynam-ics studies, achieving reliable thermodynamic propertieswhile having minimal impact on the dynamics [19]. Weuse this standard version to control the ion temperature.The electrons, however, should relax to a Fermi distri-bution, which is not possible with this classical form.We have thus produced a modified version of the Nos´e-Hoover thermostat for the dynamic electrons. It cre-ates an equilibrium distribution equal to that of a non-interacting electron gas of the same density (see Supple-mentary Information for details and derivation). With this, we ensure the ions interact with an electron sub-system with a corrected energy distribution.To demonstrate the strength of our new method, weapply it to warm dense matter (WDM). With densitiescomparable to solids and temperatures of a few electronvolts, WDM combines the need for quantum simulationsof degenerate electrons with the description of a stronglyinteracting ion component. These requirements makeWDM an ideal testbed for quantum simulations. Fur-ther, as the matter in the mantle and core of large planetsis in a WDM state [20, 21], and experiments towards in-ertial confinement fusion exhibit WDM states transientlyon the path to ignition [22, 23], simulations of WDM areof crucial importance in modern applications.Key dynamic properties of the WDM state can be rep-resented by the dynamic structure factor (DSF) [24].This quantity also connects theory and experiment:probabilities for diffraction and inelastic scattering aredirectly proportional to the DSF [25–27], allowing test-ing of WDM models. Here, we focus on the ion-ion DSFthat is defined via S ( k , ω ) = 12 πN Z exp( iωt ) h ρ ( k , t ) ρ ( − k , i dt , (2)where N is the total number of ions and ρ ( k , t ) is the spa-tial Fourier transform of the ion density. In the following,we assume the WDM system to be isotropic and spatiallyuniform. Accordingly, the structure factor depends onlyon the magnitude of the wavenumber, k = | k | . While themain contribution to S ( k, ω ) is due to direct Coulombinteractions between the ions, the modifications due toscreening are strongly affected by quantum behavior inthe electron component.Recent work has shown that predictions from standardDFT simulations for the ion-ion DSF are problematic dueto the use of the Born-Oppenheimer approximation [12].Employing a Langevin thermostat, it has been found thatthe dynamics of the electron-ion interaction may stronglychange the mode structure – in particular, the strengthof the diffusive mode. However, this approach requires avery simple, uniform frequency dependence for electron-ion collisions, which may not prove realistic in practice.It also contains an arbitrary parameter, the Langevinfriction, and is of limited predictive power as a result.We demonstrate here that our new method of Bohmiandynamics, that retains the dynamics of the electron-ioninteraction, can overcome the shortcomings of previousapproaches without introducing free parameters. Thespecific case we consider is compressed liquid aluminumwith a density of 5 . − and a temperature of 3 . k ( a ) -1 B S ( k ) ( a r b . un it s ) DFT-MDPresent Work (OFDFT) S ( k )( a . u . ) PresentWork
R=-0.013 R=0.006R=-0.005 (KS-DFT) (KS-DFT) -3 -3 -3
FIG. 2. Static ion-ion structure factors for aluminum. Thelarge graph compares our results from Bohmian dynamicswith data obtained by orbital free DFT [12] for a densityof 5 . − and a temperature of 3 . R : these values give the difference in ionic pres-sure between the methods normalized to the difference of theDFT pressures and the pressure of an ideal electron gas; thatis R = ( P Bohm − P DFT ) / ( P DFT − P ). nique described above. The comparisons with orbitalfree DFT and the computationally more intensive Kohn-Sham DFT both yield agreement within the statisticalerror of the simulations. As discussed, this match wasachieved by a single parameter fit defining λ ee . Thedifferent simulation techniques predict almost the samethermodynamics as shown by the small pressure differ-ence.Fig. 3a shows calculations of the fully frequency-dependent DSF. One can clearly notice the appearanceof side peaks in the DSF that correspond to ion acousticwaves. Their dispersion for smaller wavenumbers, andthe corresponding sound speed, are very sensitive to theinteractions present in the system. Thus, they reflect thescreening of ions by electrons as well as dynamic electron-ion collisions. For larger wavenumbers k , these modescease to exist due to increased damping. The data alsoexhibit a diffusive mode around ω = 0 although it is notas prominent as predicted in Ref. [12].The dispersion relation of the ion acoustic modes isdisplayed in Fig. 3b, which also contains results fromthe Langevin approach. The latter approach requires adhoc friction parameters which were chosen to cover therange between the classical and quantum Born limits (seeRef. [12]). The strength of the Bohmian approach lies inthe fact that it does not require a free parameter, therebyallowing access to a self-consistent prediction of the soundspeed. This comparison may also be used to assess thequality of the friction parameter applied in the Langevinapproach. For the case considered, one finds that neither a)b) k ( a ) -1 B ω (fs ) -1 ω (f s ) - FIG. 3. Results for the dynamic ion structure for aluminum at3 . . − . The upper panel shows the frequency-resolved DSF from the Bohmian dynamics. The lower panelprovides a comparison of the dispersion relation of the ionacoustic modes from our Bohmian approach with the datafrom the Langevin model of Ref. [12]. the classical nor the weak coupling Born limit are appli-cable – a finding which is typical of the WDM regime.The upper results demonstrate the strength of theBohmian approach in modeling quantum systems withstrong interactions and nonlinear ion dynamics. Forstatic and thermodynamic properties, we obtain resultsin very close agreement with DFT simulations. In ad-dition, while the standard implementation of DFT in-volves the Born-Oppenheimer approximation, our Bohmapproach can treat electrons and ions non-adiabatically,retaining the full coupling of the electron and ion dy-namics. As a result, we can investigate the changes ofthe ion modes due to dynamic electron-ion correlationswhich are inaccessible to standard DFT. In contrast toa Langevin model, we have no free parameters and canthus predict the strength of the electron drag to the ionmotion. Simulations based on time-dependent DFT [28]represent another way to avoid the Born-Oppenheimerapproximation. However, this method is numerically ex-tremely expensive, drastically limiting particle numbersand simulations times; at present, this limitation pre-cludes results for the ion modes as presented here.The principal advantage of our approach is its rela-tive numerical speed, which allows for the modeling ofquantum systems with large numbers of particles. Forcomparison, the recent time-dependent DFT simulationof Ref. [11] models a system of 128 electrons for approxi-mately 0.001 attoseconds per CPU-core and second. Thecomparative Bohmian dynamics system models 8 timesas many electrons for approximately 20 attoseconds perCPU-core and second. This vast difference in computa-tion time enables our method to access a new class ofcorrelated quantum systems. Such calculations are rele-vant not just for WDM, but also address core problemsin time-dependent chemical reactions, biological systems(e.g., protein folding), and radiation damage of materials[29–31]. Acknowledgements
This work has received support from AWE plc., theEngineering and Physical Sciences Research Council(grant numbers EP/M022331/1 and EP/N014472/1) andthe Science and Technology Facilities Council of theUnited Kingdom. This material is partially based uponwork supported by the U.S. Department of Energy, Of- fice of Science, Office of Fusion Energy Science underAward Number DE-SC0019268. c (cid:13)
British Crown Copy-right 2018/AWE.
Author contributions
This project was conceived by G.G.. The linearizationtheory, approximation scheme and numerical implemen-tation were developed by B.L. with guidance from G.G.and D.O.G.. The paper was written by B.L., D.O.G. andG.G.. Supporting calculations and theory were providedby S.R., P.M. and T.G.W..
Data availability
All data that support the findings of this study areavailable from the authors upon request.
Author declaration
The authors declare no competing financial interests. [1] Motta, M. et al.
Towards the solution of the many-electron problem in real materials: equation of stateof the hydrogen chain with state-of-the-art many-bodymethods.
Phys. Rev. X , 031059 (2017).[2] Hohenberg, P. & Kohn, W. Inhomogeneous electron gas. Phys. Rev. , B864 (1964).[3] Kohn, W. & Sham, L. J. Self-consistent equations in-cluding exchange and correlation effects.
Physical Review , A1133 (1965).[4] White, T. G., Richardson, S., Crowley, B. J. B., Pattison,L. K., Harris, J. W. O. & Gregori,G. Orbital-free density-functional theory simulations of the dynamic structurefactor of warm dense aluminum.
Phys. Rev. Lett. ,175002 (2013).[5] Knudson, M.D., Desjarlais, M.P., Becker, A., Lempke,R.W., Cochrane, K.R., Savage, M.E., Bliss, D.E., Matts-son, T.R. & Redmer, R. Direct observation of an abruptinsulator-to-metal transition in dense liquid deuterium.
Science , 1455 (2015).[6] Dharma-Wardana, M. W. C. Static and dynamicconductivity of warm dense matter within a density-functional approach: application to aluminum and gold.
Phys. Rev. E , 036401 (2006).[7] Bauernschmitt, R. & Ahlrichs, R. Treatment of electronicexcitations within the adiabatic approximation of timedependent density functional theory. Chem. Phys. Lett. , 454 (1996).[8] Burke, K. Perspective on density functional theory.
J. Chem. Phys. , 150901 (2012).[9] Ziegler, T. Approximate density functional theory asa practical tool in molecular energetics and dynamics.
Chem. Rev. , 651 (1991). [10] Jones, R. O. Density functional theory: its origins, rise toprominence, and future. Rev. Mod. Phys. , 897 (2015).[11] Baczewski, A. D., Shulenburger, L., Desjarlais, M. P.,Hansen, S. B. & Magyar, R. J. X-ray Thomson scatter-ing in warm dense matter without the Chihara decom-position. Phys. Rev. Lett. , 115004 (2016).[12] Mabey, P., Richardson, S., White, T. G., Fletcher, L. B.,Glenzer, S. H. Hartley, N. J., Vorberger, J., Gericke,D. O. & Gregori, G. A strong diffusive ion mode indense ionized matter predicted by Langevin dynamics.
Nat. Comm. , 14125 (2017).[13] Bohm, D. A suggested interpretation of the quantumtheory in terms of “hidden” variables. I. Phys. Rev. ,166 (1952).[14] Bohm, D. A collective description of electron interac-tions: III. Coulomb interactions in a degenerate electrongas. Phys. Rev. , 609 (1953).[15] Allen, M. P. & Tildesley, D. J. Computer Simulation ofLiquids. (Oxford University Press, 2017).[16] Manfredi, G. How to model quantum plasmas. arXiv:quant-ph/0505004 (2005).[17] Lado, F. Effective potential description of the quantumideal gases.
J. Chem. Phys. , 5369 (1967).[18] Dharma-Wardana, M. W. C. The classical-map hyper-netted-chain (CHNC) method and associated noveldensity-functional techniques for warm dense matter. Int. J. Quantum Chem. , 53 (2012).[19] Basconi, J. E. & Shirts, M. R. Effects of temperaturecontrol algorithms on transport properties and kineticsin molecular dynamics simulations.
J. Chem. Theory andComputation , 2887 (2013).[20] Guillot, T. Interiors of giant planets inside and outsidethe solar system. Science , 72 (1999). [21] Kerley, G. I. Equation of state and phase diagram ofdense hydrogen.
Phys. Earth and Planetary Interiors ,78 (1972).[22] Glenzer, S. H. et al. Symmetric inertial confinement fu-sion implosions at ultra-high laser energies.
Science ,1228 (2010).[23] Hu, S. X., Goncharov, V. N., Boehly, T. R., McCrory,R. L., Skupsky, S., Collins, L. A., Kress, J. D. & Militzer,B. Impact of first-principles properties of deuterium-tritium on inertial confinement fusion target designs.
Phys. Plasmas , 056304 (2015).[24] Hansen, J.P. & McDonald, I. R. Theory of Simple Liq-uids. (Academic Press, 2013).[25] Garcia Saiz, E. et al.
Probing warm dense lithium byinelastic X-ray scattering.
Nat. Physics , 940 (2008). [26] Glenzer, S. H. & Redmer, R. X-ray Thomson scatteringin high energy density plasmas. Rev. Mod. Phys. , 1625(2009).[27] Fletcher, L. B. et al. Ultrabright X-ray laser scatteringfor dynamic warm dense matter physics.
Nat. Photonics , 274 (2015).[28] Marques, M.A.L. & Gross, E.K.U. Time-dependent den-sity functional theory. Ann. Rev. Phys. Chem. , 427(2004).[29] Dill, K.A. & MacCallum, J.L. The Protein-Folding Prob-lem, 50 Years On. Science , 1042 (2012).[30] Poccia, N., Fratini, M., Ricci, A., Campi, G., Barba, L.,Vittorini-Orgeas, A., Bianconi, G., Aeppli, G. & Bian-coni, A. Evolution and control of oxygen order in acuprate superconductor.
Nat. Materials , 733 (2011).[31] Suga, M. et al. Light-induced structural changes and thesite of O=O bond formation in PSII caught by XFEL.
Nature , 131 (2017). ast Non-Adiabatic Dynamics of Many-Body Quantum Systems:Supplementary Information
B. Larder, D.O. Gericke, S. Richardson,
1, 3
P. Mabey, T. G. White, and G. Gregori ∗ Department of Physics, University of Oxford, Parks Road, Oxford OX1 3PU, UK Centre for Fusion, Space and Astrophysics, Department of Physics,University of Warwick, Coventry CV4 7AL, UK AWE, Aldermaston, Reading, Berkshire RG7 4PR, UK LULI – CNRS, Ecole Polytechnique, CEA, Universit´e Paris-Saclay, F-91128 Palaiseau Cedex, France Department of Physics, University of Nevada, Reno, Nevada 89557, USA
Bohm theory of quantum mechanics
For a single particle, the Bohm formalism follows fromwriting the wavefunction ψ as a product of amplitudeand phase information ψ ( x , t ) = R ( x , t )e iφ ( x ,t ) / ~ , (1)where x is the position vector and t the time. Using thisform for ψ , the Schr¨odinger equation yields ∂R ∂t + ∇ · (cid:18) R ∇ φm (cid:19) , (2) ∂φ∂t + V ext + V B + ( ∇ φ ) m = 0 , (3)where V ext is the external potential, and we have intro-duced the Bohm potential V B , given by V B ( x , t ) = − ~ m ∇ RR . (4)Eq. (2) provides a description of probability continuity,while Eq. (3) is of the form of a Hamilton-Jacobi equa-tion, with the additional quantum potential given by (4).Bohm’s method exploits the form of these equations bytreating the quantum system as an ensemble of classicalsystems.We can treat R as representing the density of a set ofclassical systems in configuration space. By then allow-ing the classical systems to evolve in time according tothe modified Hamilton Jacobi equation, Eq. (3), we canobtain the time evolution of R , representing the timeevolution of the probability density of the quantum sys-tem.Generalizing this idea to many-body systems, R rep-resents the probability density of classical N -body pointsin configuration space. These phase space points withinthe ensemble evolve via d x dt ≡ v = ∇ φ ( x , t ) m , (5) ∗ Corresponding author: [email protected] where x and v label the positions and velocities inconfiguration-space, respectively. This relationship tothe velocities allows for a quasi-classical description oftrajectories, or a large ensemble of N -body points to beevolved concurrently. The statistical properties of theclassical trajectories then reproduce the quantum me-chanical results. Correlation Closure
We seek a closure approximation to the full correlationfunction g appearing in Eq. (1) in the main paper, in sucha way that the Bohm potential can be constructed for thesystem. The simplest, non-trivial approximation of thiskind is a direct product representation in terms of paircorrelation functions g ( x , x , ..., x N ) ’ Y j>i g ( x i , x j ) , (6)where the pair correlation g ( x i , x j ) is the full correlationfunction with all other coordinates integrated out g ( x i , x j ) = (cid:16)Q k = i,j R d x k (cid:17) g ( x , x , ..., x N )Ω N − , (7). Here, Ω is a normalization volume.Combining Eqs. (1) of the main paper and Eq. (6) wearrive at an easily calculable expression for the Bohmpotential V B ( x | g ) ’ − ~ X j>i λ ij ( ∇ i m i + ∇ j m j ) p g ( x i , x j ) p g ( x i , x j ) . (8)For computational efficiency, we have retained interac-tions up to pair terms only, and have generalized the de-pendence on λ to a set of λ ij to accommodate differentparticle species. Fermi Statistical Corrections
At the base level of the linearized Bohmian approach,we model the system within the Hartree approximation a r X i v : . [ c ond - m a t . s t a t - m ec h ] N ov – effects due to the electron exchange interaction are notaccounted for. To rectify this, we introduce an additionalpotential term. By inverting the known pair correlationsof the non-interacting electron gas using inverse MonteCarlo (IMC) [1], we find pair potentials that reproducethe quantum correlations of the free electron gas exactly.Now we can introduce these pair potentials as additionalcontributions to the Hamiltonian.A similar procedure was pursued originally by Lado[2], and subsequently used in the classical-map method ofPerrot and Dharma-wardana for numerous applications[3–5]. This technique has also proven to be successful inrelated semi-classical correlation studies [6–8] and can beseen as analogous to correcting for the Pauli pressure inquantum hydrodynamics [9].By insisting that the total potential is equal to the IMCresult V IMC , when applied to a non-interacting electrongas, we arrive at an expression for the additional Paulipotential V P ( x ) = − V B ( x | g ) + X i,j V IMC ( x i , x j | g ) . (9)Here, g is the ideal electron gas pair correlation, and V B ( x | g ) is the Bohm potential evaluated at g . The fullsystem potential V is now given by V = V ext + V int + V B + V P , (10)where V ext is an external potential (which is set to zerothroughout our calculations), V int is the contributionfrom direct pair interactions of the simulated particles(e.g., the Coulomb potential), V B is the Bohm potential,and V P is the Pauli potential of Eq. (9). Pseudopotentials
In the description above, all particles are treated onequal footing. However, the core electrons of the ionspecies require an explicitly separate treatment from the(effectively free) valence electrons. The pseudopoten-tial approximation, commonly used in implementationsof DFT, provides a simple alternative to modeling thecore electrons directly. Core electrons are removed fromthe system to be simulated and ion-electron potentialsare constructed in a way that valence electrons form thecorrect density profile in a reference calculation. For ourpurposes, we use a reference single ion DFT calculationto determine an electron density profile and construct thepseudopotential.We use a variant of the Troullier-Martins approachalso used for DFT pseudopotentials [10]. We firstperform scalar-relativistic electron density calculationswithin DFT for an isolated ion with both n core and n core + 1 electrons. The difference of these densitiesis used as an input valence electron density. We thenconstruct the pseudopotential by requiring the followingproperties: 1. The radially integrated electron density up toa cutoff, r c , is the same for both the full ionsystem and pseudosystem. When solving forthe two-body electron-ion pseudosystem via theSchr¨odinger equation, we must have the same totalelectron density within r c as we found with the fullDFT calculation.2. The pseudopotential for r > r c is equal to the(screened) Coulomb potential3. The pseudopotential is smooth and continuous at r c (Ref. [10] discusses the specifics of the smoothingrequirement)Due to the direct connection to the approach used instandard DFT implementations, it should be possible tore-purpose existing pseudopotential generation codes forthe method explained here. This should prove usefulfor materials with complex bound state structures wheremore sophisticated pseudopotential construction schemeshave been shown to provide greater accuracy.Accompanying the application of pseudopotentials is,however, the sacrifice of correctly treating bound-freeelectron transitions. As the number of valence electronsmust be known a priori, as in implementations of orbitalfree DFT (OFDFT), we must restrict our method hereto cases where these numbers are known and reasonablywell-defined. Generalized IMC Parameter Search
We now set out how the search for the optimal valueof λ , the parameter in the linear Bohm potential, can beperformed numerically. Due to the low computationaldemands associated with calculating g for a given λ , wenote that a brute force search should be adequate in mostcases. However, a more robust approach, that is alsoapplicable to potential forms with arbitrary numbers offree parameters, is possible through a generalization ofIMC. We begin along the original lines of IMC by writingthe positional Hamiltonian of the system as a sum overpair energies: H = X α S α K α . (11)Here, K α represents the potential between two bodies ata particular distance with discrete distance bins labelledby α , and S α is the number of pairs of bodies separatedby that distance in the system, i.e., the number of pairsin the α th bin. Our goal is to obtain a set of updates toa set of parameters defining the potential, ∆ λ i , that wecan effectively adjust the mean bin counts to match ourdesired pair correlations.In our case, unlike the original formulation of IMC, thepotential K α are functionals of the thermally averagedcorrelations of the system, and are also functions of thechosen λ values, via K α = K α ( {h S α i} , { λ i } ). Here thethermal expectation h X i of an arbitrary function of thesystem coordinates X ( q ) is defined by h X i = 1 Z Z X ( q ) exp( − βH ) dq , (12)where Z is the partition function, β = 1 /k B T , and q isthe set of all positional coordinates of the system. Thatis to say, the potentials of the system depend on both theset of λ i linearization parameters, and also the thermallyaveraged set of pair correlations between bodies (as perEq. 4).Applying Eq. (12) to S α and differentiating with re-spect to λ i , we obtain ∂ h S α i ∂λ i = − β X γ M iγ ( h S α S γ i − h S α i h S γ i ) , (13)where the matrix M iγ is given by M iγ = ∂K γ ∂λ i + X δ ∂K γ ∂ h S δ i ∂ h S δ i ∂λ i . (14)Combining Eqs. (13) and (14), we find a separate set oflinear equations for each λ i M iγ = ∂K γ ∂λ i − β X (cid:15) M i(cid:15) X δ ∂K γ ∂ h S δ i ( h S δ S (cid:15) i − h S δ i h S (cid:15) i ) . (15)Using these relationships, we can now perform a directedsearch for the optimal λ i values. To leading order, we canwrite an equation for changes in h S α i in terms of changesin λ i ∆ h S α i = X i ∂ h S α i ∂λ i ∆ λ i . (16)If we perform a calculation of g with a given set of λ i values, we can simultaneously calculate values of h S α i and h S α S γ i numerically. We can then solve the linearequations (15) and, in turn, have an overdetermined setof equations (Eq. 16) for the required changes in the λ i parameters to match the desired pair correlations. Theoptimal changes in λ i can then be determined by least-squares inversion of this equation. Defining the matrixof values A αi = ∂ h S α i /∂λ i , the formal solution is∆ λ = (cid:0) A T A (cid:1) − A T ∆ h S i . (17)This can be determined directly, or implicitly throughQR factorization [11] for computational efficiency. Thisupdate procedure can then be applied iteratively to arriveat the optimal parameters λ i . Modified Thermostats
Most modern DFT-MD simulations rely on either theNos´e-Hoover or Langevin thermostat to establish an ion dynamics at a given temperature [12–14]. In the caseof the Nos´e-Hoover thermostat, an additional dynamicvariable, coupled linearly to particle momenta through afriction term, is introduced to the equations of motion[15]. The parameters for the thermostat can be chosento ensure a balance between temperature stability andthe equilibration time, but in general the dynamics areparameter-insensitive.The Langevin thermostat, on the other hand, addsboth a frictional term and a stochastic noise term tothe equations of motion [16]. Again, this ensures thecanonical distribution is sampled. The magnitude of thefriction and noise are controlled by a free parameter, theLangevin friction σ , which must be chosen to be suffi-ciently small to minimize spurious effects on the particledynamics, while remaining large enough that the correctdistribution is sampled in a reasonable time-frame. Anadvantage of this approach is that σ can be scaled upto a larger value to attempt to mimic electron-ion colli-sions, while still maintaining the canonical distribution[13]. In principle, the Born-Oppenheimer approximationcan then be employed while still attempting to includethe effects of the electron dynamics. However, the accu-rate setting of σ for this purpose is a difficult task withoutan obvious solution (Ref. [13] suggests some possible ap-proaches, although the most appropriate choice is stillunknown). Furthermore, this technique implicitly ne-glects all frequency dependence of the coupling betweenthe electron and ion dynamics as the Langevin couplingis being approximated as white noise.As our method does not rely on the Born-Oppenheimerapproximation, we can employ the parameter-insensitiveNos´e-Hoover thermostat for the ions. This represents thekey advantage of our approach: the electron dynamics ismodeled directly. Therefore, the ion thermostat – andaccordingly the ion dynamics we wish to calculate – doesnot rely on the unknown free parameter σ to mimic thedynamic electrons.In dealing directly with dynamic electron trajectories,in contrast to Born-Oppenheimer DFT-MD, we must de-velop an appropriate thermostat for the electrons as well.In the degenerate and semi-degenerate cases, the veloc-ity distribution of electrons is, however, known to be farfrom the Boltzmann distribution. While a standard ther-mostat allows us to move past the Born-Oppenheimerapproximation and include dynamic electrons, it cannotcapture the effect of Fermi-statistics on the electron mo-tion. To more accurately capture the effect of the electrondynamics on the ions, we desire a thermostat that ensuresstatic correlations are maintained, while enforcing a non-Boltzmann, that is Fermi, distribution of velocities.We highlight here modifications of standard Langevinand Nos´e-Hoover thermostats to achieve this aim. Langevin-Style Thermostat
The equations of motion of a particle in a stochasticforce field can be written quite generally as: d X = µ dt + σ · d W , (18)where X = ( X , X , X , V , V , V ) T , (19) µ = ( V , V , V , µ ( X ) , µ ( X ) , µ ( X )) T . (20)Here, µ i is the i th component of the deterministic partof the particle’s acceleration. X is a vector contain-ing the particle’s position ( X , X , X ) T and velocity( V , V , V ) T , σ represents the stochastic collision fre-quency, and W is a multidimensional Wiener process[17]. We take the driving noise to be uniform andisotropic in velocity space. Hence, σ = σ σ
00 0 σ . (21)Defining V = ( V , V , V ) T , (22) R = ( X , X , X ) T , (23) η = ( µ ( X ) , µ ( X ) , µ ( X )) T , (24)we can then write, for this system, a Fokker-Planck equa-tion for the equilibrium particle’s probability density, p ,in configuration space at long times [17]. We have ∂∂ X · ( µ p ) + σ (cid:18) ∂∂ V (cid:19) p = 0 , (25)which reduces to − V · ∂p∂ R − ∂∂ V · ( η p ) + σ (cid:18) ∂∂ V (cid:19) p = 0 . (26)We approximate the probability density distribution ofthe particles as decoupled in momentum and positionspace, via p ≡ p L = A
11 + e β ( m V / − µ c ) · e − βU ( R ) , (27)where U is the potential energy, A is a normalizationconstant, and µ c is the chemical potential. By inserting p from eq. (27) into eq. (26) we can solve for η ; theequations of motion in eq. (18) can then be written as d R = V dt (28) E ( Ha ) ρ ( E )( a . u . ) Thermostat DrivenFermi DistributionBoltzmann Distribution
FIG. 1. Reproduction of a Fermi kinetic energy distributionusing our modified thermostat. The system considered wascomprised of particles at the same density/temperature asvalence electrons in twice-compressed-Al at 3 . V = V exp( − κr ) with V = 1 Haand κ = 1 a − B was employed. Distributions are normalizedto have maxima of one. d V = F m (1 + e βE ) (cid:2) log(1 + e βE ) − βE (cid:3) dt − σ mβ (cid:18) e βE e βE (cid:19) V dt + σ · d W , (29)where m is the particle’s mass, E = m V / − µ c and F = − ∂U/∂ R . The boundary conditions were set toensure that the standard Langevin equation is recoveredat high energy. These represent the modified equationsof motion, solvable with standard methods, which havethe desired equilibrium distribution function. Nos´e-Hoover-Style Thermostat
By analogy with the Langevin case, we begin with theequations for the dynamics ˙ R = V , (30) ˙ V = F m (1 + e βE ) (cid:2) log(1 + e βE ) − βE (cid:3) − α mβ (cid:18) e βE e βE (cid:19) V , (31)where α is a new dynamic variable coupling the system tothe heat bath. In equilibrium, our deterministic systemmust now satisfy the generalized Liouville equation [18] ∂p NH ∂ R · ˙ R + ∂p NH ∂ V · ˙ V + ∂p NH ∂α ˙ α + p NH (cid:18) ∂∂ R · ˙ R + ∂∂ V · ˙ V + ∂ ˙ α∂α (cid:19) = 0 , (32)where the probability distribution p NH is given by p NH = p L exp( − βα / m α ) – in which we force α tohave a Gaussian distribution, in analogy with the stan-dard Nos´e-Hoover scheme – and m α is a thermostat mass.Again, asserting that we must recover the standard Nos´e-Hoover thermostat in the classical limit, we can solve thisequation for ˙ α . We obtain˙ αm α = − m (cid:18) e βE e βE (cid:19) (cid:20) βm V (cid:18) − e βE e βE (cid:19) + 3 (cid:21) , (33)which, in combination with Eqs. (30) and (31), formsa closed dynamical system with the desired equilibriumdistribution. For multiple particles, ˙ α is just a sum overterms of this form for each particle. Fig. 1 demonstratesthe recovery of the desired momentum distribution for asample system. Simulation Parameters
The linearization factors and Bohm potentials wereconstructed according to the theory described in the main paper. The information about static correlations neededas input, specifically the radial distribution function ofthe ionic system, was produced with OFDFT calculationrunning the using VASP package [19–22] for short times.For this purpose, we used plane wave and augmentationenergies of 250 eV and 500 eV, respectively, the 3e Alpseudopotential provided with VASP, and a simulationbox containing 256 atoms. A total of 2680 electron bandswere included, such that the occupations of the highestenergy bands were less than 10 − . This simulation wasrun for 5000 time steps of 1 fs.Electron density profiles required to produce pseudo-potentials were calculated using the GPAW DFT package[23–25] with the PBE functional [26]. The MD simulationstep was carried out with a standard MD code, inclusiveof our novel thermostat. As usual, periodic boundaryconditions were employed. The simulation box contained256 aluminum ions and the simulations were run witha time step of 0 . . u . All properties were extracted byaveraging over multiple simulation runs each having 10 time steps. [1] A.P. Lyubartsev and A. Laaksonen. Calculation of effec-tive interaction potentials from radial distribution func-tions: A reverse Monte Carlo approach, Phys. Rev. E ,3730 (1995).[2] F. Lado. Effective Potential Description of the QuantumIdeal Gases. J. Chem. Phys. , 5369 (1967).[3] M.W.C. Dharma-Wardana and F. Perrot. Simple Clas-sical Mapping of the Spin-Polarized Quantum ElectronGas: Distribution Functions and Local-Field Correc-tions, Phys. Rev. Lett. , 959 (2000).[4] M.W.C. Dharma-Wardana and M.S. Murillo. Pair-distribution functions of two-temperature two-mass sys-tems: Comparison of molecular dynamics, classical-maphypernetted chain, quantum Monte Carlo, and Kohn-Sham calculations for dense hydrogen. Phys. Rev. E ,026401 (2008).[5] M. W. C. Dharma-Wardana. The classical-map hyper-netted-chain (CHNC) method and associated noveldensity-functional techniques for warm dense matter. Int.J. Quantum Chem. , 53 (2012).[6] J. Dufty and S. Dutta. Classical representation of a quan-tum system at equilibrium: Theory.
Phys. Rev. E ,032101 (2013).[7] S. Dutta and J. Dufty. Classical representation of a quan-tum system at equilibrium: Applications. Phys. Rev. E , 032102 (2013).[8] S. Dutta and J. Dufty. Uniform electron gas at warm,dense matter conditions. Europhys. Lett. arXiv:quant-ph/0505004 (2005).[10] N. Troullier and J.L. Martins. Efficient pseudopotentialsfor plane-wave calculations. II. Operators for fast itera-tive diagonalization.
Phys. Rev. B , 8861 (1991).[11] G.H. Golub and C. F. Van Loan. Matrix Computation (Johns Hopkins University Press, 1996).[12] T.G. White, S. Richardson, B.J.B. Crowley, L. K. Patti- son, J. W O. Harris, and G. Gregori. Orbital-free density-functional theory simulations of the dynamic structurefactor of warm dense aluminum.
Phys. Rev. Lett. ,175002 (2013).[13] P. Mabey, S. Richardson, T.G. White, L. B. Fletcher,S.H. Glenzer, N. J. Hartley, J. Vorberger, D.O. Gericke,and G. Gregori. A strong diffusive ion mode in denseionized matter predicted by Langevin dynamics.
Nat.Comm. Phys. Rev. B , 144105 (2017).[15] S. Nos´e. A molecular dynamics method for simulationsin the canonical ensemble. Mol. Phys. , 255 (1984).[16] J.R. Binggeli, N. Chelikowsky. Langevin molecular dy-namics with quantum forces: Application to silicon clus-ters. Phys. Rev. B , 764 (1994).[17] G.A. Pavliotis. Stochastic Processes and Applications .(Springer, 2014).[18] G.J. Martyna and M.L. Klein. Nose-Hoover chains:thecanonical ensemble via continuous dynamics.
Phys. Rev.Lett. , 2643 (1992).[19] G. Kresse and J. Hafner. Ab initio molecular dynamicsfor liquid metals. Phys. Rev. B , 558 (1993).[20] G. Kresse and J. Hafner. Ab initio molecular-dynamicssimulation of the liquid-metalamorphous-semiconductortransition in germanium. Phys. Rev. B , 14251 (1994).[21] G. Kresse and J. Furthm¨uller. Efficient iterative schemesfor ab initio total-energy calculations using a plane-wavebasis set. Phys. Rev. B , 11169 (1996).[22] G. Kresse and J. Furthm¨uller. Efficiency of ab-initio totalenergy calculations for metals and semiconductors usinga plane-wave basis set. Comp. Mat. Science , 15 (1996).[23] J. Enkovaara, C. Rostgaard, J. J. Mortensen, J. Chen,M. Du lak, L. Ferrighi, J. Gavnholt, C. Glinsvad,V. Haikola, H. A. Hansen, H.H. Kristoffersen,M. Kuisma, A. H. Larsen, L. Lehtovaara, M. Ljung- berg, O. Lopez-Acevedo, P. G. Moses, J. Ojanen,T. Olsen, V. Petzold, N. A. Romero, J. Stausholm-Møller, M. Strange, G.A. Tritsaris, M. Vanin, M. Walter,B. Hammer, H. H¨akkinen, G. K.H. Madsen, R. M. Niem-inen, J. K. Nørskov, M. Puska, T. T. Rantala, J. Schiøtz,K. S. Thygesen, and K.W. Jacobsen. Electronic structurecalculations with GPAW: A real-space implementation ofthe projector augmented-wave method. J. Phys.: Con-densed Matter , 253202 (2010).[24] J. J. Mortensen, L. B. Hansen, and K. W. Jacobsen. Real- space grid implementation of the projector augmentedwave method. Phys. Rev. B , 035109 (2005).[25] A.H. Larsen, J.J. Mortensen, J. Blomqvist, and K.W.Jacobsen. The atomic simulation environment–a Pythonlibrary for working with atoms. J. Phys.: CondensedMatter , 273002 (2017).[26] J.P. Perdew, K. Burke, and M. Ernzerhof. GeneralizedGradient Approximation Made Simple. Phys. Rev. Lett.77