Temperature specification in atomistic molecular dynamics and its impact on simulation efficacy

Temperature is a vital thermodynamical function for physical systems. Knowledge of system temperature permits assessment of system ergodicity, entropy, system state and stability. Rapid theoretical and computational developments in the fields of condensed matter physics, chemistry, material science, molecular biology, nanotechnology and others necessitate clarity in the temperature specification. Temperature-based materials simulations, both standalone and distributed computing, are projected to grow in prominence over diverse research fields. In this article we discuss the apparent variability of temperature modeling formalisms used currently in atomistic molecular dynamics simulations, with respect to system energetics,dynamics and structural evolution. Commercial simulation programs, which by nature are heuristic, do not openly discuss this fundamental question. We address temperature specification in the context of atomistic molecular dynamics. We define a thermostat at 400K relative to a heat bath at 300K firstly using a modified ab-initio Newtonian method, and secondly using a Monte-Carlo method. The thermostatic vacancy formation and cohesion energies, equilibrium lattice constant for FCC copper is then calculated. Finally we compare and contrast the results.


Introduction
Temperature is a crucially important parameter in almost all systems and its specification is not necessarily straightforward, particularly for atomistic systems. Some simulations report temperature-based results that are within "tolerable" accuracy bounds, by relying implicitly on the temperature dependence of another parameter. For instance, some have investigated temperature effects by presuming its expansive effects on lattice parameters indirectly through empirically determined linear expansivity of the bulk material. However, motivated by the rapidly growing field of nanotechnology, it is understood that macroscopic, microscopic and nanoscopic properties are not necessarily commensurate. These emerging fields necessitate the correct thermodynamical modeling and specification of simulation temperature. The aim of this article is to compare and contrast the effects of two thermostat specifications on thermally sensitive parameters of a simple, face-centered cubic (FCC) arrangement of copper atoms. The simulations are standalone for a computationally cost-effective, small atomic array to aid comparison with respect to accuracy. Thermostats are used for a number of reasons in condensed matter simulations, such as in the search for high-temperature dynamics, annealing, formation energies, structural changes, and other effects. Additionally, limiting of steady state drifts caused by encroaching numerical errors in MD simulations, as well as provision for heat evacuation in non-equilibrium MD simulations can be done through thermostats. We examine two thermostat definitions based on

Thermostats in EPT methods
A time instant δt of simulation deterministically defines the evolution of length and energy parameters. This means that given a particular input, the application of the iterative algorithm will produce the same result at time δt. System energetics are conveyed by the Hamiltonian, where U and K are total kinetic and potential energy respectively,x is displacement and p is momentum. Instantaneous temperature T (t) is defined in terms of the instantaneous internal kinetic energy ε. Then, macroscopic system temperature equals averaged instantaneous temperature, i.e. T =⟨T (t)⟩ [1]. In the thermostatic approach variance in the Hamiltonian, as might arise from computational error accumulation, is tolerable because the thermostat is referenced to a heat bath that tightly constrains temperature variations about the equilibrium value. At equilibrium, the standard deviation σ(H) is related to isochoric heat capacity, c V , by where k B is Boltzmann constant. Internal particle motion within a system of particles is generally a complex combination of translations and rotations about a center of mass which contribute three and six degrees of freedom respectively. Computational simplifications can be achieved by subtraction of degrees of freedom when stochastic and friction are absent from the system because energy exchange between mechanically uncoupled degrees does not then occur [2,3]. If the center of mass is at rest and the particle motions are irrotational, then Using the discrete notationp i (t) :=p n i , then in the absence of friction and stochastic forces the velocity Störmer-Verlet method [4, 5, 6] discretizes particle position and velocity as Since instantaneous temperature T (t) is directly related to atomic velocities, thermostats require a mechanism to control the rate of change of particle velocities. This can be achieved either by introducing friction or by direct velocity scaling. In the first method or Langevin approach [7,8,9], Newton's second law is modified by introducing a stochastic forceR i (t) and a friction coefficient γ i (t)>0 into the acceleration of the i-th particle: In the absence of stochasticity, the friction term loses its traditional meaning and only indicates heat flow direction relative to the infinite heat bath. In which case, γ i (t)>0 implies decelerating particles and consequently loss of heat to the bath. Conversely, γ i (t)<0 implies heat gain from the bath and accelerating particles. This can also be interpreted as isochoric system pressure fluctuations. The second method or velocity scaling starts at Eq. 2 to define a scaling factor β in terms of temperatures, kinetic energies and consequently velocities [10,11,12,13] i.e. After every δt step, the velocity is therefore reassigned a valuev ′ instead ofv for the new temperature to be T ′ . Velocity scaling using constant β affects temperature distribution strongly, and a damping factor Γ∈(0,1) can be introduced [11] s.t. β = (1+Γ(T ′ /T −1)) 2 . If β is adjusted at the end of every time step δt, then dT /dt ∼ (T − T ′ ), i.e. the temperature rate of change depends on the temperature difference. The method does not remove local correlations in particle motions. Setting γ n i (t)=ξ n m i , the non-stochastic net force is where the term β=(1 − ξ n δt) describes the friction-time behavior. UsingF n i andF n+1 i in Eqs. 3 and 6 gives new position and velocity: where φ=(1 + δt 2 ξ n+1 ). Therefore, friction amounts to scaling present velocity first by β, then overall by 1/φ. For a thermostat kinetic energy is constant i.e. dε/dt=0. For the i-th particle hence the constant temperature condition: Some thermostats introduce an extra degree of freedom in the form of heat bath temperature. This conveys the extent of system-bath coupling [12,13]. Since ξ n arises from friction, the rate of system kinetic energy loss is related to the friction-time behaviour: where M >> 1 implies less frictional coupling. To determinev n+1 i in Eq.8, the term ξ n+1 must be known. Many nonlinear methods [14,15,16] arise from discretizations of Eq.12, e.g. midpoint time averaging: gives one possible thermostat when dξ/dt from Eq.12 is used. Simpler, more approximate discretizations are possible. In general, then the equations of motion and phase space can be derived. The phase space may be defined in terms of intensive macroscopic variables, which are aggregates of microscopic contributions. Explicit reference to bath temperature generally avoids systematic energy drifts, implying ergodicity [19,20,21]. In such a system, particles cannot leave the ensemble. Therefore, assuming equal Maxwell-Boltzmann particle distribution probabilities and ergodicity a priori may result in a canonical ensemble (NVT) of microstates and momenta [17]. The probability distribution of microstates over radial and momentum space takes the form Puttingp i = m iẋi allows simplification in Cartesian coordinates. The particle velocities can be shown to also follow a Maxwell-Boltzmann distribution. There are other MD thermostats that are specified under various limiting conditions [17,22,23]. Caution is needed when comparing thermostats grouped holistically within the MD approach. This is because a particular MD thermostat implementation may facilitate determination of a particular parameter, but is not generally advantageous to all situations. For instance, Hoover-Evans [20,24] and Haile-Gupta [23] thermostats are inherently time-reversible, whereas Andersen [10] and Berendsen [11] thermostats are not.

Thermostats using MC methods
A simpler thermostat can be defined non-deterministically by sampling relative particle positions using Monte-Carlo (MC) techniques [17,25]. This involves testing the effort of a random displacement of a particle within the domain (i.e. box dimensions) against the change in potential energy ∆U(x) without involving kinetic energy. The move is "acceptable" if a probability p can be found such that The microstate probability distribution ρ(x) is then This form, where particle positions are sampled exponentially without geometrical constraints, encapsulates a thermostat in the canonical ensemble that obeys Maxwell-Boltzmann (MB) thermodynamics.

Timescale and macroscopics
Timescale can be thought of as times larger than the relaxation time in interatomic particle collisions but smaller than the shortest experimental observation time for the system. The deterministic approach, unlike the MC approach, is explicit in the timescale. This naturally leads to the fundamental question of whether time-correlation functions can be found in the MC approach to describe macroscopic time-based phenomena. Instantaneous particle positions depend on the time behavior and temperature spread over particles. It is therefore important to deal with temperature-time behavior in the MC method, at the very least qualitatively, in the absence of analytical expressions [10]. The average temperature per particle is taken asT for the system in contact with a bath at temperature T bath . Many time-based system parameters of interest such as diffusion coefficients are calculated as statistical averages in the Einstein formulation, or by time-correlation integrals in the Green-Kubo formulation. At constant volume V the change in average energy is dĒ=c v dT . By Newton's law of cooling: The constant κ encapsulates temperature inhomogeneities and system shape and is estimated from the heat flow continuity and flux equations: Solving Eq. (18) givesT on a timescale such that stochasticity inT (t) averages out. From purely theoretical considerations the choice of either a specific MD, or even the MC thermostat, is subjective. However, one expects correct dynamics with a thermostat that accommodates temperature fluctuations on the timescale.

Simulation
The present work investigates the effect of the foregoing thermostat definitions, deterministic and Monte Carlo, on the equilibrium lattice constant and formation energies on a isochoric, monoatomic array of n=1583 (plus one, if with a surface adatom) copper atoms in the FCC structure, and a heat bath at 300K where necessary. The system is decidedly simple for speedier computation and comparison. The deterministic approach employs the Sutton-Chen (SC) embedded atom model (EAM) form of the Finnis-Sinclair potential [26,27,28,29]. The authors have previously reported their implementation of SC code [3], which was modified for the present work to accommodate velocity scaling-based thermostats using Störmer-Verlet integration [6], as well as the MC thermostat. EAM methods are not purely deterministic but their lower computational effort for more particles and better handling of larger timescales is attractive [30]. The EAM Hamiltonian can be written [28,30,31] where V =V (x ij ) is the pairwise interaction potential energy between particles i and j separated a distancex ij , and F =F (ρ i ) is the embedding energy of atom i as a function of the host electron density. For FCC metals the translation and rotation invariant total potential energy [4,26,32,33] can be written where and σ is a length parameter e.g. lattice constant, ε is an energy scaling factor. The fitting constants are c, m and n, where m < n. The force can be written as Hence, for the i-th particle, the force is wherex ij = (x j −x i ) is the vector between particles i and j, expressed in terms of particle positions. The Hamiltonian (H) is here accepted to within O(δt 2 ) discretization error, but is ideally constant [34]. Other temperature dependent effects like diffusion can be investigated through the effect of the thermostat definition on existing literature models. For instance, diffusion [35] is probabilistically modeled by where P m and P v are inter-lattice point migration and vacancy availability probabilities respectively, expressed in terms of migration energy E m and vacancy formation energy E v . These energies can be evaluated for each thermostat T , by considering the various vacancy formation and migratory mechanisms for the perfect crystal lattice under adatom perturbation [36]. We consider only for the {100} equivalent planes of the lattice-point to surface Schottky defect creation shown in Fig.1 and from-infinity to surface migration, shown in Fig.2. The  From-infinity to surface migration. software reported in [3] calculates the average cohesive energy, E coh , as the difference between the total potential energy of the perfect, undisturbed FCC lattice and when all N atoms are "blown apart" to infinite interatomic separation i.e.

Results and discussion
The isotherm implementations using both the deterministic model and the Monte-Carlo method were taken to be at 400K relative to a bath temperature of 300K. The two thermostats were used to calculate cohesive (E ch ), vacancy formation (E v ), Schottky formation (E sch ) (1), infinity-tosurface atom capture (E inf ) and extraction energies (E ext ), and equilibrium lattice constant (a eq ) under constant volume. Table 1 lists the input parameters for the simulation of copper. All adatom effects were simulated under the assumption that they occurred relative to planes equivalent to the (100) plane. Table 2 shows the results of the two thermostats. The point at infinity is taken approximately 100 lattice parameters away from the bulk of the array, since the atoms are assumed to be in a closed box of fixed volume.  The results shown in Table 2 were compared with literature values, where abundant empirical data exists in the literature at 400K for the deterministic approach in Method A [37,38,39,40], than for the MC approach, Method B. The MC method requires aggregation over several steps per move to give a result obeying the canonical distribution, and the actual reported results in Table 2 are time-based averages of several computations per move, with a move acceptance probability of around 30%. The move perturbations are selected randomly through the Markovchain approach, with a probability dependent on ∆x N , which represent all possible generated combinations due to perturbing inter-particle separation, x. The difficulty with the MC approach is therefore appreciable for larger ensembles owing to the total size of possible configurations. This accounts for the lower MC method accuracies depicted in the Table 2.

Conclusion
In this article we present the rationale for two broad thermostat definitions, the MD and MC thermostats. It then proceeds to describe their basic approaches and then compares and contrasts them. The choice of a thermostat for a given simulation is motivated by several factors, not least being the particular parameter of interest, rather than ease of implementation. In summary, the MD approach is generally deterministic, time-reversible, and canonical in the Hamiltonian and potential energy. The behavior of the system Hamiltonian is tolerant of drifts from the equilibrium, making the method more amenable iterative error accumulation. External degrees of freedom for the system, such as a bounding bath temperature, can readily be defined in terms of their ergodic impact on the system. The calculated dynamical trajectories can therefore be much smoother, even though method generally lacks the ability to constrain the total kinetic energy. The MC method has the main advantage of ease of thermostatic specification and description of external degrees of freedom. In this work, several difficulties were associated with the MC method, particularly that of determining time-related dynamics and the drastically increased combinations as a power of ensemble size. However, the article highlights the importance of correct temperature specification. Commercial packages are necessarily heuristic with respect to this specification and as systems evolve, particularly more towards nanoscale aggregates and ever shorter interaction timescales, the effect of temperature will become ever more important. This suggests clear directions for more future research.