Configurational thermostats for molecular systems

A new method of thermostatting non-equilibrium molecular dynamics (NEMD) simulations is described. The thermostat is based on a recently developed, entirely configurational expression for the temperature. To demonstrate this method, thermostatted NEMD simulations are performed on WCA atoms, linear, freely jointed Lennard-Jones 8-mer chains and a united-atom model of n-decane under a constant applied strain rate. The results of simulations thermostatted kinetically (the standard method) and configurationally are compared. As expected, both types of thermostat yield identical system properties for low strain rates. For higher strain rates, both thermostats yield the same qualitative dependence of system properties on applied strain rate. The great advantage of the configurational thermostat is that no a priori knowledge of the streaming velocity is required. For molecular systems and atomic systems in most flow geometries, the analytical form for the atomic streaming velocity is not known. This makes the implementation of standard kinetic thermostats highly problematic.     

Size dependence of multipolar plasmon resonance frequencies and damping rates in simple metal spherical nanoparticles

We have simulated the non-equilibrium dynamics of methanol adsorbed in FAU zeolite driven by external microwave (MW) radiation. We have modelled steady states produced by augmenting this MW-driven system with a thermostat that acts as a balancing heat sink. We have compared results from an implicit thermostat (Andersen velocity replacement) and an explicit thermostat (helium atoms subjected to Andersen velocity replacement). We find very good agreement between the implicit and explicit thermostats for energy distributions and diffusion coefficients produced under MW-heated steady-state conditions. This augurs well for the continued use of implicit thermostats, which are computationally more efficient.     

Impact of thermostat on interfacial thermal conductance prediction from non-equilibrium molecular dynamics simulations

The knowledge of interfacial thermal conductance (ITC) is key to understand thermal transport in nanostructures. The non-equilibrium molecular dynamics (NEMD) simulation is a useful tool to calculate the ITC. In this study, we investigate the impact of thermostat on the prediction of the ITC. The Langevin thermostat is found to result in larger ITC than the Nose-Hoover thermostat. In addition, the results from NEMD simulations with the Nose-Hoover thermostat exhibit strong size effect of thermal reservoirs. Detailed spectral heat flux decomposition and modal temperature calculation reveal that the acoustic phonons in hot and cold thermal reservoirs are of smaller temperature difference than optical phonons when using the Nose-Hoover thermostat, while phonons in the Langevin thermostat are of identical temperatures. Such a non-equilibrium state of phonons in the case of the Nose-Hoover thermostat reduces the heat flux of low-to-middle-frequency phonons. We also discuss how enlarging the reservoirs or adding an epitaxial rough wall to the reservoirs affects the predicted ITC, and find that these attempts could help to thermalize the phonons, but still underestimate the heat flux from low-frequency phonons.     

Energy-based coupling of smooth particle hydrodynamics and molecular dynamics with thermal fluctuations

We propose a thermodynamically consistent and energy conserving coupling scheme between the atomistic and the continuum domain. The coupling scheme links the two domains using the DPDE (Dissipative Particle Dynamics at constant Energy) thermostat and is designed to handle strong temperature gradients across the atomistic/continuum domain interface. The fundamentally different definitions of temperature in the continuum and atomistic domain – internal energy and heat capacity versus particle velocity – are accounted for in a straightforward and conceptually intuitive way by the DPDE thermostat. We verify the here proposed scheme using a fluid, which is simultaneously represented as a continuum using Smooth Particle Hydrodynamics, and as an atomistically resolved liquid using Molecular Dynamics. In the case of equilibrium contact between both domains, we show that the correct microscopic equilibrium properties of the atomistic fluid are obtained. As an example of a strong non-equilibrium situation, we consider the propagation of a steady shock-wave from the continuum domain into the atomistic domain, and show that the coupling scheme conserves both energy and shock-wave dynamics.     

不同热浴下Co团簇熔化行为的分子动力学模拟

运用分子动力学方法,采用Berendsen热浴和Nose-Hoover热浴分别研究了Co_n(n=13,55,147)团簇的熔化特性,模型采用Gupta相互作用势.模拟结果表明:两种热浴对钴团簇熔点及预熔化区间给出了基本一致的描述.所研究团簇体系在给定温度下长时间内各Co团簇中单个原子的速率(速度)分布与麦克斯韦速率(速度)分布曲线符合很好.     

不同热浴下Co团簇熔化行为的分子动力学模拟

运用分子动力学方法,采用Berendsen热浴和Nose-Hoover热浴分别研究了Con (n=13,55,147)团簇的熔化特性,模型采用Gupta相互作用势.模拟结果表明：两种热浴对钴团簇熔点及预熔化区间给出了基本一致的描述.所研究团簇体系在给定温度下长时间内各Co团簇中单个原子的速率(速度)分布与麦克斯韦速率(速度)分布曲线符合很好.     

A molecular dynamics simulation of nanoscale convective heat transfer with the effect of axial heat conduction

Effective heat dissipation from nano-fluidic devices is sometimes necessary to ensure their performance and lifespan. In the molecular dynamics simulation of nanoscale convective heat transfer, thermostats cannot be directly applied to the fluid because of the non-uniform temperature distribution. Periodic boundary is typically utilised, but unrealistic axial heat conduction exists when there is a temperature difference between the outlet and images of inlet atoms. In this paper, the effect of axial conduction caused by periodic boundary is investigated through the Péclet number (Pe). Taking viscous dissipation into consideration, the magnitude of outlet thermal diffusion is observed to decrease with increasing Pe. The local average temperature of fluid changes in an exponential form except in the region close to the outlet. Results show that the contribution of outlet axial conduction to the local average temperature is less than 2.0% when Pe > 10. The main reason is that the magnitude of fluid velocity and viscous heat dissipation in nanochannels is much larger than that in macro-channels at the same Péclet number.     

Comparing the Efficiencies of Stochastic Isothermal Molecular Dynamics Methods

Molecular dynamics typically incorporates a stochastic-dynamical device, a “thermostat,” in order to drive the system to the Gibbs (canonical) distribution at a prescribed temperature. When molecular dynamics is used to compute time-dependent properties, such as autocorrelation functions or diffusion constants, at a given temperature, there is a conflict between the need for the thermostat to perturb the time evolution of the system as little as possible and the need to establish equilibrium rapidly. In this article we define a quantity called the “efficiency” of a thermostat which relates the perturbation introduced by the thermostat to the rate of convergence of average kinetic energy to its equilibrium value. We show how to estimate this quantity analytically, carrying out the analysis for several thermostats, including the Nosé-Hoover-Langevin thermostat due to Samoletov et al. (J. Stat. Phys. 128:1321–1336, 2007) and a generalization of the “stochastic velocity rescaling” method suggested by Bussi et al. (J. Chem. Phys. 126:014101, 2007). We find efficiency improvements (proportional to the number of degrees of freedom) for the new schemes compared to Langevin Dynamics. Numerical experiments are presented which precisely confirm our theoretical estimates.     

Diffusivity and conductivity of a primitive model electrolyte in a nanopore

Equilibrium and non-equilibrium molecular dynamics simulations are applied to obtain the diffusion coefficient and electric conductivity of ions in dilute electrolytes confined in neutral cylindrical pores. The electrolyte is described with the restricted primitive model and the wall of the pore is modelled as a soft wall. The equilibrium molecular dynamics simulations show that the axial diffusion coefficient of ions decreases with increasing confinement. For a fixed pore radius the diffusion coefficient decreases with increasing number density of the ions. The current response of the system to an applied electric field is maintained at constant temperature by Gaussian isokinetic equations of motion, and at constant concentration by periodic boundary conditions with recycling of ions in the axial direction. The electric conductivity is calculated from the current density and the electric field applied for different pore sizes. In contrast to the trend in diffusivity, conductivity increases slightly in smaller pores. For a very small pore, however, conductivity is lower than the bulk, because oppositely charged ions moving in opposite directions under the electric field cannot avoid collisions with each other in a narrow channel.     

An efficient parallel algorithm for non-equilibrium molecular dynamics simulations of very large systems in planar Couette flow

In this paper we present an efficient parallel domain decomposition algorithm for non-equilibrium molecular dynamics (NEMD) simulations of large systems under planar Couette flow. We propose a modified deforming cell method that permits NEMD simulations with negligible penalties due to the Lees-Edwards periodic boundary conditions. The algorithm was used to study large systems of the Weeks-Chandler-Andersen fluid in order to obtain better viscosity results at the low shear rate regimes where the signal-to-noise ratio is very small.     

Boundary treatment effects on molecular dynamics simulations of interface thermal resistance

Molecular Dynamics simulations of heat conduction in liquid Argon confined in Silver nano-channels are performed subject to three different thermal conditions. Particularly, different surface temperatures are imposed on Silver domains using a thermostat in all and limited number of solid layers, resulting in heat flux in the liquid domain. Alternatively, energy is injected and extracted from solid layers to create a NVE liquid Argon system, which corresponds to heat flux specification. Imposition of a constant temperature region in the solid domain results in an unphysical temperature jump, indicating the presence of an artificial thermal resistance induced by the thermostat. Thermal resistance analyses for the components of each case are performed to distinguish the artificial and interface thermal resistance effects. Constant wall temperature simulations are shown to exhibit superposition of the artificial and interface thermal resistance values at the liquid/solid interface, while applying thermostat on wall layers sufficiently away from the liquid/solid interface results in consistent predictions of the interface thermal resistance. Injecting and extracting energy from each solid layer eliminates the artificial resistance. However, the method cannot directly specify a desired temperature difference between the two solid domains.     

Numerical distortion and effects of thermostat in molecular dynamics simulations of single-walled carbon nanotubes

In this paper, single-walled carbon nanotubes （SWCNTs） are studied through molecular dynamics （MD） simulation. The simulations are performed at temperatures of 1 and 300K separately, with atomic interactions characterized by the second Reactive Empirical Bond Order （REBO） potential, and temperature controlled by a certain thermostat, i.e. by separately using the velocity scaling, the Berendsen scheme, the Nose-Hoover scheme, and the generalized Langevin scheme. Results for a （5,5） SWCNT with a length of 24.5 nm show apparent distortions in nanotube configuration, which can further enter into periodic vibrations, except in simulations using the generalized Langevin thermostat, which is ascribed to periodic boundary conditions used in simulation. The periodic boundary conditions may implicitly be applied in the form of an inconsistent constraint along the axis of the nanotube. The combination of the inconsistent constraint with the cumulative errors in calculation causes the distortions of nanotubes. When the generalized Langevin thermostat is applied, inconsistently distributed errors are dispersed by the random forces, and so the distortions and vibrations disappear. This speculation is confirmed by simulation in the case without periodic boundary conditions, where no apparent distortion and vibration occur. It is also revealed that numerically induced distortions and vibrations occur only in simulation of nanotubes with a small diameter and a large length-to-diameter ratio. When MD simulation is applied to a system with a particular geometry, attention should be paid to avoiding the numerical distortion and the result infidelity.     

Nonequilibrium Statistical Mechanics and Entropy Production in a Classical Infinite System of Rotators

We analyze the dynamics of a simple but nontrivial classical Hamiltonian system of infinitely many coupled rotators. We assume that this infinite system is driven out of thermal equilibrium either because energy is injected by an external force (Case I), or because heat flows between two thermostats at different temperatures (Case II). We discuss several possible definitions of the entropy production associated with a finite or infinite region, or with a partition of the system into a finite number of pieces. We show that these definitions satisfy the expected bounds in terms of thermostat temperatures and energy flow.     

Configurational temperature thermostat for fluids undergoing shear flow: application to liquid chlorine

Non-equilibrium molecular dynamics (NEMD) simulations play a major role in characterizing the rheological properties of fluids undergoing shear flow. However, all previous studies of flows in molecular fluids either use an ‘atomic’ thermostat which makes incorrect assumptions concerning the streaming velocity of atoms within their constituent molecules, or they employ a centre of mass kinetic (COM) thermostat which only controls the temperature of relatively few degrees of freedom (3) in complex high molecular weight compounds. In the present paper we show how recently developed configurational expressions for the thermodynamic temperature can be used to develop thermostatting mechanisms which avoid both of these problems. In this work, we propose a thermostat based on a configurational expression for the temperature and apply it to NEMD simulations of chlorine undergoing Couette flow. The results so obtained are compared with those obtained using a COM kinetic thermostat. At equilibrium the properties of systems thermostatted in the two different ways are of course equivalent. We show that the two responses only differ far from equilibrium. In particular, we show that the formation of a string phase for extremely high shear rates is an artefact of the COM thermostat. At the largest shear rates studied with the configurational thermostat, no string phase is observed.     

Cell-level canonical sampling by velocity scaling for multiparticle collision dynamics simulations

A local Maxwellian thermostat for the multiparticle collision dynamics algorithm is proposed. The algorithm is based on a scaling of the relative velocities of the fluid particles within a collision cell. The scaling factor is determined from the distribution of the kinetic energy within such a cell. Thereby the algorithm ensures that the distribution of the relative velocities is given by the Maxwell–Boltzmann distribution. The algorithm is particularly useful for non-equilibrium systems, where temperature has to be controlled locally. We perform various non-equilibrium simulations for fluids in shear and pressure-driven flow, which confirm the validity of the proposed simulation scheme. In addition, we determine the dynamic structure factors for fluids with and without thermostat, which exhibit significant differences due to suppression of the diffusive part of the energy transport of the isothermal system.     

Modeling Thermostating, Entropy Currents, and Cross Effects by Dynamical Systems

A generalized multibaker map with periodic boundary conditions is shown to model boundary-driven transport, when the driving is applied by a perturbation of the dynamics localized in a macroscopically small region. In this case there are sustained density gradients in the steady state. A non-uniform stationary temperature profile can be maintained by incorporating a heat source into the dynamics, which deviates from the one of a bulk system only in a (macroscopically small) localized region such that a heat (or entropy) flux can enter an attached thermostat only in that region. For these settings the relation between the average phase-space contraction, the entropy flux to the thermostat and irreversible entropy production is clarified for stationary and non-stationary states. In addition, thermoelectric cross effects are described by a multibaker chain consisting of two parts with different transport properties, modeling a junction between two metals.     

Phonon relaxation and heat conduction in one-dimensional Fermiben Pastaben Ulam β lattices by molecular dynamics simulations

The phonon relaxation and heat conduction in one-dimensional Fermi-Pasta-Ulam (FPU) β lattices are studied by using molecular dynamics simulations. The phonon relaxation rate, which dominates the length dependence of the FPU β lattice, is first calculated from the energy autocorrelation function for different modes at various temperatures through equilibrium molecular dynamics simulations. We find that the relaxation rate as a function of wave number k is proportional to k^1.688, which leads to a N^0.41 divergence of the thermal conductivity in the framework of Green-Kubo relation. This is also in good agreement with the data obtained by non-equilibrium molecular dynamics simulations which estimate the length dependence exponent of the thermal conductivity as 0.415. Our results confirm the N^2/5 divergence in one-dimensional FPU β lattices. The effects of the heat flux on the thermal conductivity are also studied by imposing different temperature differences on the two ends of the lattices. We find that the thermal conductivity is insensitive to the heat flux under our simulation conditions. It implies that the linear response theory is applicable towards the heat conduction in one-dimensional FPU β lattices.     

A multiscale coupling method for the modeling of dynamics of solids with application to brittle cracks

We present a multiscale model for numerical simulations of dynamics of crystalline solids. The method combines the continuum nonlinear elasto-dynamics model, which models the stress waves and physical loading conditions, and molecular dynamics model, which provides the nonlinear constitutive relation and resolves the atomic structures near local defects. The coupling of the two models is achieved based on a general framework for multiscale modeling – the heterogeneous multiscale method (HMM). We derive an explicit coupling condition at the atomistic/continuum interface. Application to the dynamics of brittle cracks under various loading conditions is presented as test examples.     

Computing the local pressure in molecular dynamics simulations

Computer simulations of inhomogeneous soft matter systems often require accurate methods for computing the local pressure. We present a simple derivation, based on the virial relation, of two equivalent expressions for the local (atomistic) pressure in a molecular dynamics simulation. One of these expressions, previously derived by other authors via a different route, involves summation over interactions between particles within the region of interest; the other involves summation over interactions across the boundary of the region of interest. We illustrate our derivation using simulations of a simple osmotic system; both expressions produce accurate results even when the region of interest over which the pressure is measured is very small.