NonEquilibriumHeatExchange

class NonEquilibriumHeatExchange(configuration, heating_power, heat_source, heat_sink, exchange_interval=1, update_profile_interval=0, profile_resolution=None)

A class that implements a heat flow by constant heating power exchange technique via a hook function.

Parameters:
  • configuration (BulkConfiguration) – The initial configuration on which the heat flow simulation will be performed.
  • heating_power (PhysicalQuantity of type energy/time) – The thermal energy that is transferred from the cold to the hot region per time.
  • heat_source (str | list of ints) – The tag or list of indices defining the group of atoms in the hot region.
  • heat_sink (str | list of ints) – The tag or list of indices defining the group of atoms in the cold region.
  • exchange_interval (int) – The interval to perform energy exchange by rescaling the velocities.
    Default: 1 (every step).
  • update_profile_interval (int) – The interval at which a measurement of the temperature profile is performed and added to the average profile. Set to zero to disable on-the-fly measurement.
    Default: 0
  • profile_resolution (PhysicalQuantity of type length) – The spatial resolution for the on-the-fly calculation of the temperature profile.
    Default: 2.0*Angstrom
callInterval()
Returns:The call interval of this hook function.
Return type:int
heatCurrent()
Returns:The average heat current.
Return type:PhysicalQuantity of type energy/time
nlprint(stream=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='UTF-8'>)

Print a formatted string with the average temperature profile, as well as the average heat current.

Parameters:stream (A stream that supports strings being written to using 'write'.) – The stream the temperature profile is written to.
temperatureProfile()
Returns:The temperature profile as a tuple containing bin-centers and corresponding temperature values.
Return type:(PhysicalQuantity of type length, PhysicalQuantity of type temperature)

Usage Example

Perform a non-equilibrium molecular dynamics simulation with a constant heating power on a silicon crystal with a single germanium layer, to obtain the thermal conductance through this impurity layer.

# Add tags
bulk_configuration.addTags('heat_sink',   [
    108, 109, 110, 111, 112, 113, 114, 115, 116, 225, 226, 227, 228, 229, 230, 231, 232, 233, 342,
    343, 344, 345, 346, 347, 348, 349, 350, 459, 460, 461, 462, 463, 464, 465, 466, 467, 576, 577,
    578, 579, 580, 581, 582, 583, 584, 693, 694, 695, 696, 697, 698, 699, 700, 701, 810, 811, 812,
    813, 814, 815, 816, 817, 818, 927, 928, 929, 930, 931, 932, 933, 934, 935,
])
bulk_configuration.addTags('heat_source', [
    0,   1,   2,   3,   4,   5,   6,   7,   8, 117, 118, 119, 120, 121, 122, 123, 124, 125, 234,
    235, 236, 237, 238, 239, 240, 241, 242, 351, 352, 353, 354, 355, 356, 357, 358, 359, 468, 469,
    470, 471, 472, 473, 474, 475, 476, 585, 586, 587, 588, 589, 590, 591, 592, 593, 702, 703, 704,
    705, 706, 707, 708, 709, 710, 819, 820, 821, 822, 823, 824, 825, 826, 827,
])

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------

potentialSet = Tersoff_SiGe_1989()
calculator = TremoloXCalculator(parameters=potentialSet)
bulk_configuration.setCalculator(calculator)

# -------------------------------------------------------------
# Molecular Dynamics
# -------------------------------------------------------------

initial_velocity = MaxwellBoltzmannDistribution(
    temperature=300.0*Kelvin,
    remove_center_of_mass_momentum=True
)

method = NVEVelocityVerlet(
    time_step=1*femtoSecond,
    initial_velocity=initial_velocity,
)

# Use a heating power of 50 nW.
heating_power = 50.0e-9*Joule/Second
momentum_exchange_hook = NonEquilibriumHeatExchange(
    configuration=bulk_configuration,
    heating_power=heating_power,
    heat_source='heat_source',
    heat_sink='heat_sink',
    update_profile_interval=10,
    profile_resolution=2.0*Ang
)

md_trajectory = MolecularDynamics(
    bulk_configuration,
    constraints=[],
    trajectory_filename='Si_w_Ge_layer_NEMD.hdf5',
    steps=100000,
    log_interval=500,
    post_step_hook=momentum_exchange_hook,
    method=method
)

# Print the information from the NEMD heat flow simulation.
nlprint(momentum_exchange_hook)

Here, only the molecular dynamics block of the script is shown. The full script can be found found in the file nonequilibriumheatexchange.py.

Notes

  • The NonEquilibriumHeatExchange method can be used to run MolecularDynamics simulations with a constant heat flux between the two selected reservoirs. The technique can therefore be used to simulate the thermal conductance of a given configuration (see e.g, [JJ99]).

  • The reverse non-equilibrium heat exchange method is implemented as a class, which can be used as a post_step_hook in MolecularDynamics.

  • The velocities in the two reservoir regions are rescaled in such a way that a well-defined amount kinetic energy is added to the atoms in the heat_source region whereas the same amount is taken away from the kinetic energy in the heat_sink region. The magnitude of the kinetic energy that is transferred per time from the heat_sink to the heat_source can be specified by the heating_power parameter. By default the rescaling is performed at every MD step, if desired, the interval can be given by the user via the parameter exchange_interval.

  • If a non-zero update_profile_interval is specified, a temperature profile is stored on-the-fly during the simulation. After the MD run, the profiles may be accessed using the method temperatureProfile.

  • The thermal bulk conductivity can be obtained by plotting the temperature profile of the system using the MD-Analyzer tool and fitting its gradient \(\nabla T\) along the transport direction. The thermal conductivity is then obtained via Fourier’s law

    \[\kappa = - \frac{\dot{Q}}{A \nabla T}\]

where \(\dot{Q}\) is the specified heating_power and A is the cross-sectional area perpendicular to the transport direction.

  • The thermal conductance across an interface G (Kapitza conductance) can be obtained as

    \[G = -\frac{\dot{Q}}{\Delta T}\]

where \(\Delta T\) is the temperature drop across the interface.

[JJ99]Philippe Jund and Rémi Jullien. Molecular-dynamics calculation of the thermal conductivity of vitreous silica. Phys. Rev. B, 59:13707–13711, 1999. doi:10.1103/PhysRevB.59.13707.