Author Topic: How to Apply Heat Flux in NEMD Simulations Using the "Exchange Interval" Option  (Read 76784 times)

0 Members and 1 Guest are viewing this topic.

Offline Priyadarshini

  • Regular QuantumATK user
  • **
  • Posts: 7
  • Country: in
  • Reputation: 0
    • View Profile
Hi all,

I'm setting up a non-equilibrium molecular dynamics (NEMD) simulation in QuantumATK to calculate the thermal conductivity of a material. I understand that QuantumATK allows heat flux to be imposed by exchanging kinetic energy between two regions (hot and cold) using the "Exchange Interval" parameter.

However, I'm unclear about a few implementation details:

How exactly is the heat flux imposed using the exchange method in QuantumATK?
Is the energy exchanged directly between predefined regions at every exchange interval step? How does this correspond to a physical heat flux value?
How should I define the size and location of the heat source and sink regions?
Should they be single layers of atoms, or should they span a certain number of unit cells?
How do I correctly choose the exchange interval?
What criteria should I consider when selecting the interval value to ensure physical accuracy and numerical stability?
How can I extract or calculate the actual heat flux value that is being applied during the simulation?
Is there a way to monitor the energy exchanged over time to derive the heat flux?

Any guidance would be greatly appreciated!

Offline Anders Blom

  • QuantumATK Staff
  • Supreme QuantumATK Wizard
  • *****
  • Posts: 5743
  • Country: dk
  • Reputation: 112
    • View Profile
    • QuantumATK at Synopsys
The heat sinks need to be a bit bigger than just a few atoms, usually 4-5 layers is probably fine.
The other fundamental questions should be answered in the references under https://docs.quantumatk.com/manual/Types/NonEquilibriumMomentumExchange/NonEquilibriumMomentumExchange.html
« Last Edit: September 3, 2025, 03:17 by Anders Blom »

Offline Julian Schneider

  • QuantumATK Staff
  • QuantumATK Guru
  • *****
  • Posts: 166
  • Country: dk
  • Reputation: 25
    • View Profile
The RNMED method is implemented as described in these papers:

F. Mueller-Plathe. A simple nonequilibrium molecular dynamics method for calculating the thermal conductivity. J. Chem. Phys., 106(14):6082–6085, 1997https://doi.org/10.1063/1.473271.
C. Nieto-Draghi and J. B. Avalos. Non-equilibrium momentum exchange algorithm for molecular dynamics simulation of heat flow in multicomponent systems. Mol. Phys., 101(14):2303–2307, 2003.https://doi.org/10.1080/0026897031000154338
In short, the energy is exchanged by swapping the momenta of the hottest atom in the heat sink region and the coldest atom in the heat source region. This way an external heat flow between sink and source region is established which is counterbalanced by a heat flux in the configuration from source to sink. When steady-state is reached the external heat flux (which can be measured from the exchanged energy) is equal to the internal heat flux. The thermal conductivity can then be calculated from this heat flux and the temperature gradient using Fourier's law, as described in the tutorial.

We normally recommend a slab-like setup where the top and bottom layers of the surfaces are frozen (after equilibration) and thin heat sink and source regions are assigned adjacent to these frozen caps, e.g. as shown in the attached picture.
These regions can be specified as tags. The thickness should be between 5 and 10 Ang.
The exchange interval should be chosen as big as possible to avoid a too large temperature gradient, and but not too big so that one still gets a measurable gradient. Recommended values are 100 -200. If the temperature gradient becomes too large, then one should increase the value. I have attached an example script of a thermal conductivity simulation.
The heat flux can be extracted when opening the MD-Trajectory in the MD-Analyzer. The measurement "average_heat_current" can be enabled which show the time evolution of the average energy exchanged per time. This can be used to monitor if one has reached a steady-state. The heat flux is obtained by dividing the heat current by the cross-sectional area of the configuration, i.e. the A-B-surface area if the current is along C. The temperature gradient can be obtained from the MD-Analyzer, by plotting the temperature profile and fitting the slope in a selected region as shown in the third attached picture. The tutorial https://docs.quantumatk.com/tutorials/interfacial_thermal_conductance/interfacial_thermal_conductance.html#investigating-the-grain-size-dependence-of-the-thermal-conductivity outlines this, but the functionality displayed there is a bit outdated.

Since this gives you full control over all details, but can be a bit cumbersome to do, we have implemented a convenience function that does combines all these steps and prints the thermal conductivity to the log: https://docs.quantumatk.com/manual/Types/ThermalConductivity/ThermalConductivity.html#thermalconductivity

Offline Priyadarshini

  • Regular QuantumATK user
  • **
  • Posts: 7
  • Country: in
  • Reputation: 0
    • View Profile
Hello QuantumATK team,

Thank you for your guidance; the process is now clear to me. I go through the reference image attached showing a temperature profile analysis for SiO₂ quartz (MTP pretrained, RNEMD) using the Molecular Dynamics Analyzer in QuantumATK.
Could you please share the official documentation or example workflow/script that explains how to set up and perform the RNEMD (Reverse Non-Equilibrium Molecular Dynamics) simulation for SiO₂ using the MTP potential and generate and fit the temperature profile as shown in this plot?

Any references or sample Python scripts for automating this process in the QuantumATK Script Generator would be very helpful.

Thank you in advance for your guidance!