NPTMartynaTobiasKlein¶

Included in QATK.Dynamics

class NPTMartynaTobiasKlein(initial_velocity=None, time_step=None, reservoir_temperature=None, thermostat_timescale=None, reservoir_pressure=None, barostat_timescale=None, coupling_mask=None, heating_rate=None, compression_rate=None, chain_length=None)¶

The NPTMartynaTobiasKlein integrator class which implements the barostat proposed by Martyna, Tobias, and Klein.

Parameters:

initial_velocity (ConfigurationVelocities | ZeroVelocities | MaxwellBoltzmannDistribution) – A class that implements a distribution of initial velocities for the particles in the MD simulation.
Default: MaxwellBoltzmannDistribution
time_step (PhysicalQuantity of type time) – The time-step interval used in the MD simulation.
Default: 1 * fs
reservoir_temperature (PhysicalQuantity of type temperature | list) – The reservoir temperature in the simulation.
Default: 300.0 * Kelvin
thermostat_timescale (PhysicalQuantity of type time) – The time constant for MTK temperature coupling.
Default: 100.0 * fs
reservoir_pressure (PhysicalQuantity of type pressure) – The reservoir pressure in the simulation. The pressure can be given either as a scalar value, as a vector of length 3 denoting the hydrostatic components, as a vector of length 6 in Voigt notation, or as a 3x3 tensor. A scalar value will result in isotropic pressure coupling, whereas for all other representations, each pressure component will be coupled to the barostat independently.
Default: 1.0 * bar
barostat_timescale (PhysicalQuantity of type time) – The time constant for MTK pressure coupling.
Default: 1000.0 * fs
coupling_mask (array of bools) – The mask determining which elements of the stress tensor, the barostat should couple to. Can be given as a vector of length 3, denoting the hydrostatic components, as a vector of length 6 in Voigt notation, or as a 3x3 tensor. By default the barostat couples only to the diagonal elements.
Default: True for diagonal elements, False for shear components.
heating_rate (PhysicalQuantity of type temperature/time | None) – The heating rate of the target temperature. A value of None disables the heating of the system.
Default: None
compression_rate (PhysicalQuantity of type pressure/time | None) – The compression rate of the target pressure.
Default: 0.0 * bar / fs
chain_length (int) – The number of subsequent thermostats.
Default: 3

static getLogInfo(step, time, md_quantities)¶

Get the log information for the MD step.

Parameters:

step (int) – The current MD step.
time (PhysicalQuantity of type time) – The current MD time in fs.
md_quantities (MDQuantities) – The MD quantities to log.

Returns:

The header rows and the row with the MD information.

Return type:

list of str, str

isotropicCoupling()¶

Returns:: Whether the barostat applies isotropic or anisotropic pressure coupling.
Return type:: bool

kineticEnergy(configuration)¶

Parameters:: configuration (DistributedConfiguration) – The current configuration to calculate the kinetic energy of.
Returns:: The kinetic energy of the current configuration.
Return type:: PhysicalQuantity of type energy

kineticEnergyLeapfrog(configuration, forces)¶

Parameters:

configuration (DistributedConfiguration) – The current configuration to calculate the kinetic energy of.
forces (PhysicalQuantity of type energy per distance) – The current forces on the configuration.

Returns:

The kinetic energy of the current configuration.

Return type:

PhysicalQuantity of type energy

thermostats()¶

Returns:: The list of thermostats.
Return type:: list

timeStep()¶

Returns:: The time step.
Return type:: PhysicalQuantity of type time

uniqueString()¶: Return a unique string representing the state of the object.

Usage Example¶

Perform a molecular dynamics run of 50 steps on FCC Si, using the MartynaTobiasKlein barostat with isotropic pressure coupling:

# Set up configuration
bulk_configuration = BulkConfiguration(
    bravais_lattice=FaceCenteredCubic(5.4306*Angstrom),
    elements=[Silicon, Silicon],
    fractional_coordinates=[[0.0, 0.0, 0.0], [0.25, 0.25, 0.25]]
    )

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

# Set up a MTK-barostat with isotropic pressure coupling.
method = NPTMartynaTobiasKlein(
    time_step=1*femtoSecond,
    reservoir_temperature=300*Kelvin,
    thermostat_timescale=100*femtoSecond,
    reservoir_pressure=1.0*bar,
    barostat_timescale=500.0*femtoSecond,
    heating_rate=0*Kelvin/picoSecond,
    compression_rate=0.0*bar/picoSecond,
    chain_length=3,
    initial_velocity=None
)

# Run MD simulation
md_trajectory = MolecularDynamics(
    bulk_configuration,
    constraints=[],
    trajectory_filename='trajectory.nc',
    steps=50,
    log_interval=10,
    method=method
)

Apply anisotropic pressure coupling only to the zz-direction of the pressure tensor.

method = NPTMartynaTobiasKlein(
    time_step=1*femtoSecond,
    reservoir_temperature=300*Kelvin,
    thermostat_timescale=100*femtoSecond,
    reservoir_pressure=[0.0, 0.0, 10.0]*GPa,
    barostat_timescale=500.0*femtoSecond,
    heating_rate=0*Kelvin/picoSecond,
    coupling_mask=[False, False, True, False, False, False],
    compression_rate=0.0*bar/picoSecond,
    chain_length=3,
    initial_velocity=None
)

Notes¶

The MartynaTobiasKlein barostat can be used to simulate a canonical NPT-ensemble. Similar to the NVTNoseHoover algorithm, chains of subsequent thermostats and barostats are used to suppress large temperature and pressure oscillations. The equations of motion are taken from Ref. [1], whereas the integration algorithm is implemented as described in Ref [2].
If a scalar reservoir_pressure value is given, isotropic pressure coupling is applied, i.e. all cell vectors are rescaled by the same factor, and the isotropic pressure \(P=1/3 (P_{xx} + P_{yy} + P_{zz})\) is used.

If a pressure vector (of length 3 or 6, in Voigt notation), or a 3x3-tensor is given as reservoir_pressure, anisotropic pressure coupling is applied, meaning that each degree of freedom of the cell vectors will be rescaled independently.

You can enable/disable the coupling to selected components of the pressure tensor via the parameter coupling_matrix. This parameter only has an effect with anisotropic pressure coupling.

By default, the barostat couples only to the hydrostatic pressure, i.e. the diagonal elements of the pressure tensor.

If a non-zero heating_rate is specified, the reservoir temperature will be changed linearly during the simulation, according to the specified heating rate.

If a non-zero compression_rate is specified, the reservoir pressure will be changed linearly during the simulation, according to the specified compression rate.
This method is only supported in the general dynamics frameworks. See the documentation on the SoftMatterDynamicsSimulation for details on which methods are supported in each framework.