NPTBerendsen

class NPTBerendsen(initial_velocity=None, time_step=None, reservoir_temperature=None, thermostat_timescale=None, reservoir_pressure=None, barostat_timescale=None, compressibility=None, coupling_mask=None, heating_rate=None, compression_rate=None)

NPT Berendsen integrator class which implements the Berendsen barostat.

Parameters:
  • initial_velocity (ConfigurationVelocities | ZeroVelocities | MaxwellBoltzmannDistribution) – A class that implements a distribution of initial velocities for the particles in the MD simulation.
    Default: ZeroVelocities
  • time_step (PhysicalQuantity of type time) – The time-step interval used in the MD simulation.
    Default: 1.0 * fs
  • reservoir_temperature (PhysicalQuantity of type temperature | list) – The reservoir temperature in the simulation. The temperature can be given as a single temperature value for the entire system, or as a list of 2-tuples of str and PhysicalQuantity of type temperature, applying local thermostats to the tagged groups of atoms. E.g. [('group1', 280.0 * Kelvin), ('group2', 320.0 * Kelvin)].
    Default: 300.0 * Kelvin
  • thermostat_timescale (PhysicalQuantity of type time) – The time constant for Berendsen 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 Berendsen pressure coupling.
    Default: 1000.0 * fs
  • compressibility (PhysicalQuantity of type pressure**-1) – The estimated compressibility of the system relating volume changes to pressures changes.
    Default: 1.0e-4*bar**-1
  • 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
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
thermostats()
Returns:The list of thermostats.
Return type:list
timeStep()
Returns:The time step.
Return type:PhysicalQuantity of type time

Usage Example

Perform a molecular dynamics run of 50 steps on FCC Si, using the Berendsen 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 Berendsen barostat with isotropic pressure coupling.
method = NPTBerendsen(
    time_step=1*femtoSecond,
    reservoir_temperature=300*Kelvin,
    thermostat_timescale=100*femtoSecond,
    reservoir_pressure=1.0*bar,
    barostat_timescale=500.0*femtoSecond,
    compressibility=1.0e-4*bar**-1,
    heating_rate=0*Kelvin/picoSecond,
    compression_rate=0.0*bar/picoSecond,
    initial_velocity=None
)

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

nptberendsen.py

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

method = NPTBerendsen(
    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,
    compressibility=1.0e-4*bar**-1,
    heating_rate=0*Kelvin/picoSecond,
    coupling_mask=[False, False, True, False, False, False],
    compression_rate=0.0*bar/picoSecond,
    initial_velocity=None
)

nptberendsen2.py

Notes

  • The Berendsen barostat can be used to approximately simulate a canonical NPT-ensemble, as described in Ref. [BPvG+84].

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

  • The compressibility parameter is used to determine the relationship between volume and pressure changes, when the simulation cell is rescaled, which affects how well the observed barostat time scale agrees with the specified barostat_timescale parameter. However, if not set to extreme and unphysical values, the compressibility parameter should not influence the result of the simulation.

[BPvG+84]H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak. Molecular dynamics with coupling to an external bath. J. Chem. Phys., 81(8):3684–3690, 1984. doi:http://dx.doi.org/10.1063/1.448118.