NonEquilibriumMomentumExchange

class NonEquilibriumMomentumExchange(configuration, exchange_interval, heat_source, heat_sink, update_profile_interval=0, profile_resolution=None)

A class that implements a heat flow by non-equilibrium momentum exchange (reverse non-equilibrium MD (RNEMD)) technique via a hook function.

Parameters:
  • configuration (BulkConfiguration) – The initial configuration on which the heat flow simulation will be performed.
  • exchange_interval (int) – After each exchange_interval iterations, a velocity swapping is performed between the hot and the cold region.
  • heat_source (str | list of ints | tuple of type (float, PhysicalQuantity of type length)) – The heat source atoms can be defined by giving a tag that specifies which atoms are in the heat source, a list of atomic indices, or by defining a spatial region. The spatial region is defined as 2-tuple of the origin and length. The origin is the fractional position along the c-axis and the length is the thickness along the c-axis.
  • heat_sink (str | list of ints | tuple of type (float, PhysicalQuantity of type length)) – The heat sink atoms can be defined by giving a tag that specifies which atoms are in the heat sink, a list of atomic indices, or by defining a spatial region. The spatial region is defined as 2-tuple of the origin and length. The origin is the fractional position along the c-axis and the length is the thickness along the c-axis.
  • 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
averageHeatCurrent()
Returns:The average heat current.
Return type:PhysicalQuantity of type energy/time
callInterval()
Returns:The call interval of this hook function.
Return type:int
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 molecular dynamics run of 5000 steps on a carbon nanotube. Momenta of two particles in the hot, respectively cold region are exchanged every 100 steps.

# Set up configuration
bulk_configuration = NanoTube(3,1)

# Add tags
nano_tube.addTags('cold', [0, 1, 2, 3, 5, 6, 7, 10, 11, 16, 17])
nano_tube.addTags('hot',  [24, 30, 31, 35, 36, 37, 41, 42, 45, 46, 47, 48, 49, 50, 51])

# Set calculator
calculator = TremoloXCalculator(parameters=Tersoff_C_2010())
bulk_configuration.setCalculator(calculator)

method = NVEVelocityVerlet(
    time_step=0.5*femtoSecond,
    initial_velocity=MaxwellBoltzmannDistribution()
)

momentum_exchange_hook = NonEquilibriumMomentumExchange(
    configuration=bulk_configuration,
    exchange_interval=100,
    hot_region='hot',
    cold_region='cold',
    update_profile_interval=0,
    resolution=2.0*Ang
)

md_trajectory = MolecularDynamics(
    bulk_configuration,
    constraints=[],
    trajectory_filename='trajectory.nc',
    steps=5000,
    log_interval=50,
    post_step_hook=momentum_exchange_hook,
    method=method
)

nonequilibriummomentumexchange.py

Perform a non-equilibrium momentum exchange simulation on a polyvinyl chloride polymer system. The heat source atoms are determined at each exchange step as any atoms within 10 Angstrom of the origin of the c-axis and the heat sink atoms are any atoms within 10 Angstroms of the middle of the cell.

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

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

momentum_exchange_hook = NonEquilibriumMomentumExchange(
    configuration=bulk_configuration,
    exchange_interval=100,
    heat_source=(0.0, 10.0*Angstrom),
    heat_sink=(0.5, 10.0*Angstrom),
    update_profile_interval=0,
    profile_resolution=2.0*Ang
)

constraints = [FixCenterOfMass()]

md_trajectory = MolecularDynamics(
    bulk_configuration,
    constraints=constraints,
    trajectory_filename='nemd_polymer.hdf5',
    steps=50000,
    log_interval=100,
    measurement_hook=momentum_exchange_hook,
    method=method
)

bulk_configuration = md_trajectory.lastImage()
nlprint(momentum_exchange_hook)

nemd_polymer.py

Notes

  • The reverse non-equilibrium momentum-exchange (RNEMD) method is implemented as a class, which can be used as a measurement_hook in MolecularDynamics.

  • The average heat current is recorded as a MD measurement. The time trace of this value can be used to monitor the convergence of the calculation. It can be easily accessed by opening the MD trajectory using the Movie Tool in the Builder. It may also be accessed by using the measurement method of MDTrajectory. For example:

    times, average_heat_currents = md_trajectory.measurement('average_heat_current')
    
  • The heat sink and heat source atoms can be selected in several different ways. The first way to is explicitly specify the indices of the atoms. The second way is to specify the names of tagged groups of atoms. The third way is to define a region of space whose atoms will be included. The third option assumes that the heat current is measured along the c-axis of the simulation cell and is defined in terms of an origin and a length. The origin is given in fractional coordinates and the length is given as a PhysicalQuantity with length units. This method is ideal for working with liquid systems where atoms may drift from their initial positions. See the descriptions of the arguments heat_source and heat_sink above for details.

  • The momenta of the coldest particle in the heat_source group of atoms and the hottest particle in the heat_sink group are exchanged at a given interval. This exchange produces a temperature gradient, and an internal flow of energy between the two groups of atoms. Both quantities can be measured and the thermal conductance can be calculated using Fourier’s law. For more details regarding the RNEMD technique, please consult Refs. [MP97] and [NDA03].

  • 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.

[MP97]F. Mueller-Plathe. A simple nonequilibrium molecular dynamics method for calculating the thermal conductivity. J. Chem. Phys., 106(14):6082–6085, 1997. doi:10.1063/1.473271.
[NDA03]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. doi:10.1080/0026897031000154338.