Non-equilibrium momentum exchange algorithm for molecular dynamics simulation of heat flow in multicomponent systems

Green-Kubo expressions for transport coefficients from dissipative particle dynamics simulations revisited

This article addresses the debate about the correct application of Green-Kubo expressions for transport coefficients from dissipative particle dynamics simulations. We demonstrate that the Green-Kubo expressions are valid provided that (i) the dynamic model conserves the physical property, whose transport is studied, and (ii) the fluctuations satisfy detailed balance. As a result, the traditional expressions used in molecular dynamics can also be applied to dissipative particle dynamics simulations. However, taking the calculation of the shear viscosity as a paradigmatic example, a random contribution, whose strength scales as 1/δt1/2, with δt the time-step, can cause difficulties if the stress tensor is not separated into the different contributions. We compare our expression to that of Ernst and Brito (M. H. Ernst and R. Brito, Europhys. Lett., 2006, 73, 183-189), which arises from a diametrically different perspective. We demonstrate that the two expressions are completely equivalent and find exactly the same result both analytically and numerically. We show that the differences are not due to the lack of time-reversibility but instead from a pre-averaging of the random contributions. Despite the overall validity of Green-Kubo expressions, we find that the Einstein-Helfand relations (D. C. Malaspina et al. Phys. Chem. Chem. Phys., 2023, 25, 12025-12040) do not suffer from the need to decompose the stress tensor and can readily be used with a high degree of accuracy. Consequently, Einstein-Helfand relations should be seen as the preferred method to calculate transport coefficients from dissipative particle dynamics simulations.

Predicting thermodiffusion in simple binary fluid mixtures

The predictive capabilities of some existing theoretical models to quantify thermodiffusion have been investigated in this work. To do so, the tests have been performed on two model fluids, the hard-sphere and the Lennard-Jones (including spheres and dimers) ones, exploring different mixtures and thermodynamic conditions thanks to extensive molecular simulations. It has been confirmed that the thermal diffusion factor should be expressed as the sum of one term related to the isotope effect and one term related to the "chemical" effects and that a kinetic term is required to quantify thermodiffusion from the gas state to the liquid state. In addition, regarding the isotope effects, it has been obtained that none of the available theoretical models are able to yield a reasonable prediction relatively to the molecular simulations results and that the moment of inertia contribution is one order of magnitude smaller than the mass contribution in the liquid state. Finally, concerning the chemical effects, it has been shown the Shukla and Firoozabadi model, complemented with a kinetic term, is probably the most reasonable option to estimate the chemical contribution to the thermal diffusion factor, even if it fails in capturing the effect of the asymmetry in size and in shape between the species. Overall, this works confirms that there is still a lack of a generic model able to predict accurately thermal diffusion factors, or equivalently Soret coefficient, in simple binary mixtures from the gas state to the liquid state.

Thermal properties of single-layer MoS2-WS2 alloys enabled by machine-learned interatomic potentials

Two-dimensional (2D) quantum materials are poised to transform conventional electronics for a wide spectrum of applications that will encompass chemical sciences. For the study of thermal transport in single-layer (1L)...

In silico all-atom approach to thermodiffusion in dilute aqueous solutions

Thermodiffusion (or thermophoresis) is the phenomenon by which the spatial distributions of constituents of liquid or gas phases become inhomogeneous in response to a temperature gradient. It has been evidenced in a variety of systems and has many practical applications as well as implications in the context of the origins of life. A complete molecular picture of thermophoresis is still missing, and phenomenological approaches are often employed to account for the experimental observations. In particular, the amplitude of the resulting concentration-gradients (quantified by the Soret coefficient) depends on many factors that are not straightforwardly rationalized. All-atom molecular dynamics simulations appear as an exquisite tool to shed light on the molecular origins for this phenomenon in molecular systems, but the practical implementation of thermophoretic settings in silico poses significant challenges. Here, we propose a robust approach to tackle thermophoresis in dilute realistic solutions at the molecular level. We rely on a recent enhanced heat-exchange algorithm to generate temperature-gradients. We carefully assess the convergence of thermophoretic simulations in dilute aqueous solutions. We show that simulations typically need to be propagated on long timescales (hundreds of nanoseconds). We find that the magnitude of the temperature gradient and the box sizes have little effect on the measured Soret coefficients. Practical guidelines are derived from such observations. Provided with this reliable setup, we discuss the results of thermophoretic simulations on several examples of molecular, neutral solutes, which we find in very good agreement with experimental measurements regarding the concentration-, mass-, and temperature-dependence of the Soret coefficient.

Hierarchically hydrogen-bonded graphene/polymer interfaces with drastically enhanced interfacial thermal conductance

Interfacial thermal transport is a critical physical process determining the performance of many material systems with small-scale features. Recently, self-assembled monolayers and polymer brushes have been widely used to engineer material interfaces presenting unprecedented properties. Here, we demonstrate that poly(vinyl alcohol) (PVA) monolayers with hierarchically arranged hydrogen bonds drastically enhance interfacial thermal conductance by a factor of 6.22 across the interface between graphene and poly(methyl methacrylate) (PMMA). The enhancement is tunable by varying the number of grafted chains and the density of hydrogen bonds in the unique hierarchical hydrogen bond network. The extraordinary enhancement results from a synergy of hydrogen bonds and other structural and thermal factors including molecular morphology, chain orientation, interfacial vibrational coupling and heat exchange. Two types of hydrogen bonds, i.e. PVA-PMMA hydrogen bonds and PVA-PVA hydrogen bonds, are analyzed and their effects on various structural and thermal properties are systematically investigated. These results are expected to provide new physical insights for interface engineering to achieve tunable thermal management and energy efficiency in a wide variety of systems involving polymers and biomaterials.

Chemical assembling of amine functionalized boron nitride nanotubes onto polymeric nanofiber film for improving their thermal conductivity

The interfacial effect between an organic matrix and inorganic nanofillers on the thermal conductivity of a polymer composite was systematically explored by assembling amine-functionalized boron nitride nanotubes onto an electrospun nanofiber.

Polymeric Self-Assembled Monolayers Anomalously Improve Thermal Transport across Graphene/Polymer Interfaces

Ultralow thermal conductivities of bulk polymers greatly limit their applications in areas demanding fast heat dissipation, such as flexible electronics and microelectronics. Therefore, polymeric composites incorporating highly thermally conductive filler materials (e.g., graphene and carbon nanotubes) have been produced to address the issue. However, despite some enhancement, thermal conductivities of these materials are still far below theoretical predictions, mainly due to the inefficient thermal transport across material interfaces. Here, using molecular dynamics simulations, we demonstrate that polyethylene (PE) self-assembled monolayer (SAM) functionalized graphene surfaces at a high grafting density can drastically improve interfacial thermal conduction between graphene and the matrix of poly(methyl methacrylate) (PMMA). In contrast to abrupt temperature drop across pristine graphene/PMMA interfaces, temperature field in the vicinity of a PE-grafted graphene/PMMA interface is continuous with a smoother transition and higher thermal conductance. This anomalous improvement is attributed to three factors that closely relate to the grafting density of the SAM of PE. First, the SAM with high grafting densities features highly extended chains that enhance along-chain thermal conduction. Second, the strong covalent bonding between the SAM and the graphene facilitates heat transfer at their joints. Third, the SAM and the PMMA matrix are both organic materials, leading to enhanced interfacial vibrational coupling. Molecular mechanisms underpinning these phenomena are systematically elucidated by analyzing the temperature field, density distribution, Herman's orientation factor, the vibrational density of states, cumulative correlation factor, the integrated autocorrelation of interfacial heat power, and interfacial adhesion. All results suggest the incorporation of SAMs at a high grafting density or extremely extended polymer brushes for drastically improved interfacial thermal transport between hard and soft materials toward a wide range of applications.

Study of thermal conductivity of ice clusters after impact deposition on the silica surfaces using the ReaxFF reactive force field

During aircraft or spacecraft missions, ice accumulates on different parts of their surface elements. An important parameter affecting the ability to remove this ice from the surface is the heat transfer characteristics of the accumulated ice. The ice heat transfer is related to the process of ice formation and its density and internal structure. In this study we investigate the effects of the ice and silica structure and the ice cluster attachment mechanism to the silica surface on the thermal conductivity (TC) of the attached ice cluster using the ReaxFF reactive force field. The purpose of this study is to investigate the thermal transport in amorphous and crystalline ice after high-velocity deposition on the silica surfaces. A dual thermostat method has been applied for the calculation of TC values. The validity of this method has been verified by comparing the calculated values of TC for crystal and amorphous ice with available experimental values. Our calculations show that the TC values of both crystal and amorphous ice drop after deposition on the silica surfaces. This decrease in the TC is more significant for the ice deposition on suboxide silica surfaces. Furthermore, crystal ice shows higher TC values than amorphous ice after accumulation. However, when crystal ice impacts on the silica surface at 1 km s(-1) impact speed, the crystalline shape of the ice cluster is lost to a considerable level and the TC values obtained for the ice clusters in such cases are closer to amorphous ice TC values. We observed a decrease in the TC values when ionic species are added inside the ice clusters.

Tuning thermal conductivity of crystalline polymer nanofibers by interchain hydrogen bonding

Interchain hydrogen bonds enhance thermal conduction in crystalline polymer nanofibers by confining torsional motion of polymer chains and by increasing the group velocity of phonons.

Effects of the amino acid sequence on thermal conduction through β-sheet crystals of natural silk protein

Molecular mechanisms underpinning the thermal transport process through three types of β-sheets are studied to reveal the intrinsic sequence effects.

Hydrogen bonding-assisted thermal conduction in β-sheet crystals of spider silk protein

Using atomistic simulations, we demonstrate that β-sheet, an essential component of spider silk protein, has a thermal conductivity 1-2 orders of magnitude higher than that of some other protein structures reported in the literature. In contrast to several other nanostructured materials of similar bundled/layered structures (e.g. few-layer graphene and bundled carbon nanotubes), the β-sheet is found to uniquely feature enhanced thermal conductivity with an increased number of constituting units, i.e. β-strands. Phonon analysis identifies inter-β-strand hydrogen bonding as the main contributor to the intriguing phenomenon, which prominently influences the state of phonons in both low- and high-frequency regimes. A thermal resistance model further verifies the critical role of hydrogen bonding in thermal conduction through β-sheet structures.

Molecular dynamics algorithm enforcing energy conservation for microcanonical simulations

A reversible algorithm [enforced energy conservation (EEC)] that enforces total energy conservation for microcanonical simulations is presented. The key point is the introduction of the discrete-gradient method to define the forces from the conservative potentials, instead of the direct use of the force field at the actual position of the particle. We have studied the performance and accuracy of the EEC in two cases, namely Lennard-Jones fluid and a simple electrolyte model. Truncated potentials that usually induce inaccuracies in energy conservation are used. In particular, the reaction field approach is used in the latter. The EEC is able to preserve energy conservation for a long time, and, in addition, it performs better than the Verlet algorithm for these kinds of simulations.

Si/Ge superlattice nanowires with ultralow thermal conductivity

The engineering of nanostructured materials with very low thermal conductivity is a necessary step toward the realization of efficient thermoelectric devices. We report here the main results of an investigation with nonequilibrium molecular dynamics simulations on thermal transport in Si/Ge superlattice nanowires aiming at taking advantage of the inherent one dimensionality and the combined presence of surface and interfacial phonon scattering to yield ultralow values for their thermal conductivity. Our calculations revealed that the thermal conductivity of a Si/Ge superlattice nanowire varies nonmonotonically with both the Si/Ge lattice periodic length and the nanowire cross-sectional width. The optimal periodic length corresponds to an order of magnitude (92%) decrease in thermal conductivity at room temperature, compared to pristine single-crystalline Si nanowires. We also identified two competing mechanisms governing the thermal transport in superlattice nanowires, responsible for this nonmonotonic behavior: interface modulation in the longitudinal direction significantly depressing the phonon group velocities and hindering heat conduction, and coherent phonons occurring at extremely short periodic lengths counteracting the interface effect and facilitating thermal transport. Our results show trends for superlattice nanowire design for efficient thermoelectrics.

A low-frequency wave motion mechanism enables efficient energy transport in carbon nanotubes at high heat fluxes

The great majority of investigations of thermal transport in carbon nanotubes (CNTs) in the open literature focus on low heat fluxes, that is, in the regime of validity of the Fourier heat conduction law. In this paper, by performing nonequilibrium molecular dynamics simulations we investigated thermal transport in a single-walled CNT bridging two Si slabs under constant high heat flux. An anomalous wave-like kinetic energy profile was observed, and a previously unexplored, wave-dominated energy transport mechanism is identified for high heat fluxes in CNTs, originated from excited low frequency transverse acoustic waves. The transported energy, in terms of a one-dimensional low frequency mechanical wave, is quantified as a function of the total heat flux applied and is compared to the energy transported by traditional Fourier heat conduction. The results show that the low frequency wave actually overtakes traditional Fourier heat conduction and efficiently transports the energy at high heat flux. Our findings reveal an important new mechanism for high heat flux energy transport in low-dimensional nanostructures, such as one-dimensional (1-D) nanotubes and nanowires, which could be very relevant to high heat flux dissipation such as in micro/nanoelectronics applications.

Molecular dynamics simulations of thermal transport in porous nanotube network structures

Carbon nanotube based 3D nanostructures have shown a lot of promise towards designing next generation of multi-functional systems, such as nano-electronic devices. Motivated by their recent successful experimental synthesis as well as characterization, and realizing that thermal dissipation is an important concern in proposed devices because of ever-increasing power density, we have investigated the phononic thermal transport behavior in 3D porous nanotube network structures using reverse non-equilibrium molecular dynamics simulations. Based on our study, the length scale associated with the distance between nanotube junctions emerges as the most dominating parameter that governs phonon scattering (hence the characteristic mean free path) and the heat flow in these nanostructures at molecular length scales. However, because of their spatial inhomogeneity, we show that the aerial density of carbon nanotubes (normal to heat flow) is also of critical importance in determining their system-level thermal conductivity. Based on our findings, we postulate that both parameters should be considered while designing nano-devices where thermal management is relevant.

Thermal conductivity of ionic systems from equilibrium molecular dynamics

Thermal conductivities of ionic compounds (NaCl, MgO, Mg(2)SiO(4)) are calculated from equilibrium molecular dynamics simulations using the Green-Kubo method. Transferable interaction potentials including many-body polarization effects are employed. Various physical conditions (solid and liquid states, high temperatures, high pressures) relevant to the study of the heat transport in the Earth's mantle are investigated, for which experimental measures are very challenging. By introducing a frequency-dependent thermal conductivity, we show that important coupled thermoelectric effects occur in the energy conduction mechanism in the case of liquid systems.

A uniform source-and-sink scheme for calculating thermal conductivity by nonequilibrium molecular dynamics

A uniform source-and-sink (USS) scheme, which combines features of the reverse [F. Müller-Plathe, J. Chem. Phys. 106, 6082 (1997)] and improved relaxation [B. Y. Cao, J. Chem. Phys. 129, 074106 (2008)] methods, is developed to calculate the thermal conductivity by nonequilibrium molecular dynamics (NEMD). The uniform internal heat source and sink are realized by exchanging the velocity vectors of individual atoms in the right half and left half systems, and produce a periodically quadratic temperature profile throughout the system. The thermal conductivity can be easily extracted from the mean temperatures of the right and left half systems rather than by fitting the temperature profiles. In particular, this scheme greatly increases the relaxation of the exited localized phonon modes which often worsen the calculation accuracy and efficiency in most other NEMD methods. The calculation of the thermal conductivities of solid argon shows that the simple USS scheme gives accurate results with fast convergence.

Thermal conductivity of periclase (MgO) from first principles

We combine first-principles calculations of forces with the direct nonequilibrium molecular dynamics method to determine the lattice thermal conductivity k of periclase (MgO) up to conditions representative of the Earth's core-mantle boundary (136 GPa, 4100 K). We predict the logarithmic density derivative a=(∂ln k/∂ln ρ)(T)=4.6±1.2 and that k=20±5  Wm(-1) K(-1) at the core-mantle boundary, while also finding good agreement with extant experimental data at much lower pressures.

Influence of the interaction potential and of the temperature on the thermodiffusion (Soret) coefficient in a model system

In this paper the thermodiffusive behavior of an equimolar binary mixture subject to repulsive potentials of the form (sigma/r)(n) is investigated by using nonequilibrium molecular dynamics (NEMD) and the thermodiffusion (Soret) coefficient, S(T), is computed in a wide range of temperatures. With the aim to contribute to the study of the dependence of the Soret coefficient on the interaction potential, the exponent n of the potential is varied from 1 to 12, that is from a pseudocoulombian to a pseudohard-sphere interaction. The steady state equation is integrated for the composition function under reasonable assumptions and it is shown that in some cases the request for it to be linear cannot be satisfied. For this reason nonlinear functions are used to fit the NEMD composition data. The simulations indicate a negligible dependence of S(T) on the composition (in the composition range here considered) while the dependence on the temperature is more marked. The computed values of S(T) as a function of the temperature are fitted with analytical functions. It is found that with n> or =3 (medium and short range interaction) the model system behaves like a dilute gas mixture with the Soret coefficient varying with the temperature almost like 1/T. In the case of n=1 (long range interaction), S(T) has a more complex dependence on T: in particular it shows a change of sign. The analytical fitting functions, S(T)(T), are used in the integrated steady state equation thus obtaining the steady state composition profile and its comparison with the NEMD results indicates the grounding of the approach here proposed.

System size and control parameter effects in reverse perturbation nonequilibrium molecular dynamics

The issue of system size effects in the reverse perturbation nonequilibrium molecular dynamics method for determining transport coefficients of fluids is examined for the case of the Lennard-Jones model. It is found that when adequate precautions are observed in obtaining linear temperature or momentum profiles, a 250 atom system is adequate for determining the thermal conductivity and the shear viscosity. Also, a means of determining the uncertainties in the transport coefficients is described. The conclusion is that this method is computationally competitive with other simulation methods for estimating transport coefficients.

Transport coefficients and dynamic properties of hydrogen sulfide from molecular simulation

Molecular-dynamics simulation results on thermodynamic and transport properties of pure H2S under conditions of practical interest are presented. Our data are in very good quantitative agreement with the scarce experimental data and estimates on thermophysical properties of this substance. Our results serve as a test of the validity of the intermolecular potential used in the simulations as well as the consistency of the existing data in the studied range. New simulation data on thermal conductivity at low temperature as well as in supercritical states are also reported. Furthermore, we present a comparative analysis between the local order in the liquid phase of pure hydrogen sulfide and water, due to the molecular analogies between both substances, and its relation with the formation of H=S bonds. Our results indicate that under the same corresponding thermodynamic states, H2S is a much less structured substance, with a first solvation shell with a dodecahedral order instead of the tetrahedral order observed in water.

Thermodynamic and transport properties of carbon dioxide from molecular simulation

Monte Carlo and molecular dynamics simulations have been used in order to test the ability of a three center intermolecular potential for carbon dioxide to reproduce literature experimental thermophysical values. In particular, both the shear viscosity under supercritical conditions and along the phase coexistence line, as well as the thermal conductivity under supercritical conditions, have been calculated. Together with the already reported excellent agreement for the phase coexistence densities, the authors find that the agreement with experimental values is, in general, good, except for the thermal conductivity at low density. Although extended versions of the model were employed, which include an explicit account of bending and vibrational degrees of freedom, a significant difference was still found with respect to the reported experimental value.

Dynamical and structural properties of benzene in supercritical water

We have employed an anisotropic united atom model of benzene (R. O. Contreras, Ph.D. thesis, Universitat Rovira i Virgili 2002) that reproduces the quadrupolar moment of this molecule through the inclusion of seven point charges. We show that this kind of interaction is required to reproduce the solvation of these molecules in supercritical water. We have computed self-diffusion coefficient and Maxwell-Stefan coefficients as well as the shear viscosity for the mixture water-benzene at supercritical conditions. A strong density and composition dependence of these properties is observed. In addition, our simulations are in qualitative agreement with the experimental evidence that, at medium densities (0.6 g/cm(3) and 673 K), almost half of the benzene molecules have one hydrogen bond with water molecules. We also observe that these bonds are longer lived than the corresponding hydrogen bonds between water molecules. Similarly, we obtain an important reduction of the dielectric constant of the mixture with the increment of the amount of benzene molecules at medium and high densities.

Thermal Conductivities of Molecular Liquids by Reverse Nonequilibrium Molecular Dynamics

The reverse nonequilibrium molecular dynamics method for thermal conductivities is adapted to the investigation of molecular fluids. The method generates a heat flux through the system by suitably exchanging velocities of particles located in different regions. From the resulting temperature gradient, the thermal conductivity is then calculated. Different variants of the algorithm and their combinations with other system parameters are tested: exchange of atomic velocities versus exchange of molecular center-of-mass velocities, different exchange frequencies, molecular models with bond constraints versus models with flexible bonds, united-atom versus all-atom models, and presence versus absence of a thermostat. To help establish the range of applicability, the algorithm is tested on different models of benzene, cyclohexane, water, and n-hexane. We find that the algorithm is robust and that the calculated thermal conductivities are insensitive to variations in its control parameters. The force field, in contrast, has a major influence on the value of the thermal conductivity. While calculated and experimental thermal conductivities fall into the same order of magnitude, in most cases the calculated values are systematically larger. United-atom force fields seem to do better than all-atom force fields, possibly because they remove high-frequency degrees of freedom from the simulation, which, in nature, are quantum-mechanical oscillators in their ground state and do not contribute to heat conduction.

Computing the Soret coefficient in aqueous mixtures using boundary driven nonequilibrium molecular dynamics

We have computed the Soret coefficient in aqueous mixtures using a boundary driven nonequilibrium molecular dynamics algorithm and standard molecular force fields. The choice of this specific approach is justified by the nature of the mixtures studied here. Four aqueous solutions, including methanol, ethanol, acetone, and dimethyl-sulfoxide (DMSO) have been studied at ambient conditions for different compositions. The experimental behavior of water-alcohol mixtures was reproduced, including the change of sign of the Soret coefficient with composition, in excellent agreement with existing experimental data. The methodology has been applied to obtain pure predictions for water-acetone and water-DMSO where no experimental data are accessible. A change of sign is also observed in the same range of composition as in water-alcohol mixtures. It is suggested that the nature and strength of the molecular interactions, rather than the mass or shape ratio of the components, dominates the behavior of the Soret coefficient versus composition for the aqueous associating mixtures studied here.