md.integrate¶
Overview
md.integrate.berendsen |
Applies the Berendsen thermostat. |
md.integrate.brownian |
Brownian dynamics. |
md.integrate.langevin |
Langevin dynamics. |
md.integrate.mode_standard |
Enables a variety of standard integration methods. |
md.integrate.mode_minimize_fire |
Energy Minimizer (FIRE). |
md.integrate.npt |
NPT Integration via MTK barostat-thermostat. |
md.integrate.nph |
NPH Integration via MTK barostat-thermostat.. |
md.integrate.nve |
NVE Integration via Velocity-Verlet |
md.integrate.nvt |
NVT Integration via the Nosé-Hoover thermostat. |
Details
Integration methods.
To integrate the system forward in time, an integration mode must be set. Only one integration mode can be active at
a time, and the last integrate.mode_*
command before the hoomd.run()
command is the one that will take effect. It is
possible to set one mode, run for a certain number of steps and then switch to another mode before the next run
command.
The most commonly used mode is :py:class`mode_standard`. It specifies a standard mode where, at each time
step, all of the specified forces are evaluated and used in moving the system forward to the next step.
:py:class`mode_standard` doesn’t integrate any particles by itself, one or more compatible integration methods must
be specified before the staring a hoomd.run()
. Like commands that specify forces, integration methods are
persistent and remain set until they are disabled.
To clarify, the following series of commands will run for 1000 time steps in the NVT ensemble and then switch to NVE for another 1000 steps:
all = group.all()
integrate.mode_standard(dt=0.005)
nvt = integrate.nvt(group=all, kT=1.2, tau=0.5)
run(1000)
nvt.disable()
integrate.nve(group=all)
run(1000)
You can change integrator parameters between runs:
integrator = integrate.nvt(group=all, kT=1.2, tau=0.5)
run(100)
integrator.set_params(kT=1.0)
run(100)
This code snippet runs the first 100 time steps with kT=1.2 and the next 100 with kT=1.0.
-
class
hoomd.md.integrate.
berendsen
(group, kT, tau)¶ Applies the Berendsen thermostat.
Parameters: - group (
hoomd.group
) – Group to which the Berendsen thermostat will be applied. - kT (
hoomd.variant
orfloat
) – Temperature of thermostat. (in energy units). - tau (float) – Time constant of thermostat. (in time units)
berendsen
rescales the velocities of all particles on each time step. The rescaling is performed so that the difference in the current temperature from the set point decays exponentially: Berendsen et. al. 1984.\[\frac{dT_\mathrm{cur}}{dt} = \frac{T - T_\mathrm{cur}}{\tau}\]Attention
berendsen
does not function with MPI parallel simulations.Attention
berendsen
does not integrate rotational degrees of freedom.-
disable
()¶ Disables the integration method.
Examples:
method.disable()
Executing the disable command will remove the integration method from the simulation. Any
hoomd.run()
command executed after disabling an integration method will not apply the integration method to the particles during the simulation. A disabled integration method can be re-enabled withenable()
.
-
randomize_velocities
(seed)¶ Assign random velocities and angular momenta to particles in the group, sampling from the Maxwell-Boltzmann distribution. This method considers the dimensionality of the system and particle anisotropy, and removes drift (the center of mass velocity).
New in version 2.3.
Parameters: seed (int) – Random number seed Note
Randomization is applied at the start of the next call to
hoomd.run()
.Example:
integrator = md.integrate.berendsen(group=group.all(), kT=1.0, tau=0.5) integrator.randomize_velocities(seed=42) run(100)
- group (
-
class
hoomd.md.integrate.
brownian
(group, kT, seed, dscale=False, noiseless_t=False, noiseless_r=False)¶ Brownian dynamics.
Parameters: - group (
hoomd.group
) – Group of particles to apply this method to. - kT (
hoomd.variant
orfloat
) – Temperature of the simulation (in energy units). - seed (int) – Random seed to use for generating \(\vec{F}_\mathrm{R}\).
- dscale (bool) – Control \(\lambda\) options. If 0 or False, use \(\gamma\) values set per type. If non-zero, \(\gamma = \lambda d_i\).
- noiseless_t (bool) – If set true, there will be no translational noise (random force)
- noiseless_r (bool) – If set true, there will be no rotational noise (random torque)
brownian
integrates particles forward in time according to the overdamped Langevin equations of motion, sometimes called Brownian dynamics, or the diffusive limit.\[ \begin{align}\begin{aligned}\frac{d\vec{x}}{dt} = \frac{\vec{F}_\mathrm{C} + \vec{F}_\mathrm{R}}{\gamma}\\\langle \vec{F}_\mathrm{R} \rangle = 0\\\langle |\vec{F}_\mathrm{R}|^2 \rangle = 2 d k T \gamma / \delta t\\\langle \vec{v}(t) \rangle = 0\\\langle |\vec{v}(t)|^2 \rangle = d k T / m\end{aligned}\end{align} \]where \(\vec{F}_\mathrm{C}\) is the force on the particle from all potentials and constraint forces, \(\gamma\) is the drag coefficient, \(\vec{F}_\mathrm{R}\) is a uniform random force, \(\vec{v}\) is the particle’s velocity, and \(d\) is the dimensionality of the system. The magnitude of the random force is chosen via the fluctuation-dissipation theorem to be consistent with the specified drag and temperature, \(T\). When \(kT=0\), the random force \(\vec{F}_\mathrm{R}=0\).
brownian
generates random numbers by hashing together the particle tag, user seed, and current time step index. See C. L. Phillips et. al. 2011 for more information.Attention
Change the seed if you reset the simulation time step to 0. If you keep the same seed, the simulation will continue with the same sequence of random numbers used previously and may cause unphysical correlations.
For MPI runs: all ranks other than 0 ignore the seed input and use the value of rank 0.
brownian
uses the integrator from I. Snook, The Langevin and Generalised Langevin Approach to the Dynamics of Atomic, Polymeric and Colloidal Systems, 2007, section 6.2.5, with the exception that \(\vec{F}_\mathrm{R}\) is drawn from a uniform random number distribution.In Brownian dynamics, particle velocities are completely decoupled from positions. At each time step,
brownian
draws a new velocity distribution consistent with the current set temperature so thathoomd.compute.thermo
will report appropriate temperatures and pressures if logged or needed by other commands.Brownian dynamics neglects the acceleration term in the Langevin equation. This assumption is valid when overdamped: \(\frac{m}{\gamma} \ll \delta t\). Use
langevin
if your system is not overdamped.You can specify \(\gamma\) in two ways:
- Use
set_gamma()
to specify it directly, with independent values for each particle type in the system. - Specify \(\lambda\) which scales the particle diameter to \(\gamma = \lambda d_i\). The units of \(\lambda\) are mass / distance / time.
brownian
must be used with integrate.mode_standard.kT can be a variant type, allowing for temperature ramps in simulation runs.
A
hoomd.compute.thermo
is automatically created and associated with group.Examples:
all = group.all(); integrator = integrate.brownian(group=all, kT=1.0, seed=5) integrator = integrate.brownian(group=all, kT=1.0, dscale=1.5) typeA = group.type('A'); integrator = integrate.brownian(group=typeA, kT=hoomd.variant.linear_interp([(0, 4.0), (1e6, 1.0)]), seed=10)
-
disable
()¶ Disables the integration method.
Examples:
method.disable()
Executing the disable command will remove the integration method from the simulation. Any
hoomd.run()
command executed after disabling an integration method will not apply the integration method to the particles during the simulation. A disabled integration method can be re-enabled withenable()
.
-
set_gamma
(a, gamma)¶ Set gamma for a particle type.
Parameters: set_gamma()
sets the coefficient \(\gamma\) for a single particle type, identified by name. The default is 1.0 if not specified for a type.It is not an error to specify gammas for particle types that do not exist in the simulation. This can be useful in defining a single simulation script for many different types of particles even when some simulations only include a subset.
Examples:
bd.set_gamma('A', gamma=2.0)
-
set_gamma_r
(a, gamma_r)¶ Set gamma_r for a particle type.
Parameters: set_gamma_r()
sets the coefficient \(\gamma_r\) for a single particle type, identified by name. The default is 1.0 if not specified for a type. It must be positive or zero, if set zero, it will have no rotational damping or random torque, but still with updates from normal net torque.Examples:
bd.set_gamma_r('A', gamma_r=2.0) bd.set_gamma_r('A', gamma_r=(1,2,3))
-
set_params
(kT=None)¶ Change langevin integrator parameters.
Parameters: kT ( hoomd.variant
orfloat
) – New temperature (if set) (in energy units).Examples:
integrator.set_params(kT=2.0)
- group (
-
class
hoomd.md.integrate.
langevin
(group, kT, seed, dscale=False, tally=False, noiseless_t=False, noiseless_r=False)¶ Langevin dynamics.
Parameters: - group (
hoomd.group
) – Group of particles to apply this method to. - kT (
hoomd.variant
orfloat
) – Temperature of the simulation (in energy units). - seed (int) – Random seed to use for generating \(\vec{F}_\mathrm{R}\).
- dscale (bool) – Control \(\lambda\) options. If 0 or False, use \(\gamma\) values set per type. If non-zero, \(\gamma = \lambda d_i\).
- tally (bool) – (optional) If true, the energy exchange between the thermal reservoir and the particles is
tracked. Total energy conservation can then be monitored by adding
langevin_reservoir_energy_groupname
to the logged quantities. - noiseless_t (bool) – If set true, there will be no translational noise (random force)
- noiseless_r (bool) – If set true, there will be no rotational noise (random torque)
Translational degrees of freedom
langevin
integrates particles forward in time according to the Langevin equations of motion:\[ \begin{align}\begin{aligned}m \frac{d\vec{v}}{dt} = \vec{F}_\mathrm{C} - \gamma \cdot \vec{v} + \vec{F}_\mathrm{R}\\\langle \vec{F}_\mathrm{R} \rangle = 0\\\langle |\vec{F}_\mathrm{R}|^2 \rangle = 2 d kT \gamma / \delta t\end{aligned}\end{align} \]where \(\vec{F}_\mathrm{C}\) is the force on the particle from all potentials and constraint forces, \(\gamma\) is the drag coefficient, \(\vec{v}\) is the particle’s velocity, \(\vec{F}_\mathrm{R}\) is a uniform random force, and \(d\) is the dimensionality of the system (2 or 3). The magnitude of the random force is chosen via the fluctuation-dissipation theorem to be consistent with the specified drag and temperature, \(T\). When \(kT=0\), the random force \(\vec{F}_\mathrm{R}=0\).
langevin
generates random numbers by hashing together the particle tag, user seed, and current time step index. See C. L. Phillips et. al. 2011 for more information.Attention
Change the seed if you reset the simulation time step to 0. If you keep the same seed, the simulation will continue with the same sequence of random numbers used previously and may cause unphysical correlations.
For MPI runs: all ranks other than 0 ignore the seed input and use the value of rank 0.
Langevin dynamics includes the acceleration term in the Langevin equation and is useful for gently thermalizing systems using a small gamma. This assumption is valid when underdamped: \(\frac{m}{\gamma} \gg \delta t\). Use
brownian
if your system is not underdamped.langevin
uses the same integrator asnve
with the additional force term \(- \gamma \cdot \vec{v} + \vec{F}_\mathrm{R}\). The random force \(\vec{F}_\mathrm{R}\) is drawn from a uniform random number distribution.You can specify \(\gamma\) in two ways:
- Use
set_gamma()
to specify it directly, with independent values for each particle type in the system. - Specify \(\lambda\) which scales the particle diameter to \(\gamma = \lambda d_i\). The units of \(\lambda\) are mass / distance / time.
langevin
must be used withmode_standard
.kT can be a variant type, allowing for temperature ramps in simulation runs.
A
hoomd.compute.thermo
is automatically created and associated with group.Warning
When restarting a simulation, the energy of the reservoir will be reset to zero.
Examples:
all = group.all(); integrator = integrate.langevin(group=all, kT=1.0, seed=5) integrator = integrate.langevin(group=all, kT=1.0, dscale=1.5, tally=True) typeA = group.type('A'); integrator = integrate.langevin(group=typeA, kT=hoomd.variant.linear_interp([(0, 4.0), (1e6, 1.0)]), seed=10)
-
disable
()¶ Disables the integration method.
Examples:
method.disable()
Executing the disable command will remove the integration method from the simulation. Any
hoomd.run()
command executed after disabling an integration method will not apply the integration method to the particles during the simulation. A disabled integration method can be re-enabled withenable()
.
-
set_gamma
(a, gamma)¶ Set gamma for a particle type.
Parameters: set_gamma()
sets the coefficient \(\gamma\) for a single particle type, identified by name. The default is 1.0 if not specified for a type.It is not an error to specify gammas for particle types that do not exist in the simulation. This can be useful in defining a single simulation script for many different types of particles even when some simulations only include a subset.
Examples:
bd.set_gamma('A', gamma=2.0)
-
set_gamma_r
(a, gamma_r)¶ Set gamma_r for a particle type.
Parameters: set_gamma_r()
sets the coefficient \(\gamma_r\) for a single particle type, identified by name. The default is 1.0 if not specified for a type. It must be positive or zero, if set zero, it will have no rotational damping or random torque, but still with updates from normal net torque.Examples:
langevin.set_gamma_r('A', gamma_r=2.0) langevin.set_gamma_r('A', gamma_r=(1.0,2.0,3.0))
-
set_params
(kT=None, tally=None)¶ Change langevin integrator parameters.
Parameters: - kT (
hoomd.variant
orfloat
) – New temperature (if set) (in energy units). - tally (bool) – (optional) If true, the energy exchange between the thermal reservoir and the particles is
tracked. Total energy conservation can then be monitored by adding
langevin_reservoir_energy_groupname
to the logged quantities.
Examples:
integrator.set_params(kT=2.0) integrator.set_params(tally=False)
- kT (
- group (
-
class
hoomd.md.integrate.
mode_minimize_fire
(dt, Nmin=5, finc=1.1, fdec=0.5, alpha_start=0.1, falpha=0.99, ftol=0.1, wtol=0.1, Etol=1e-05, min_steps=10, group=None, aniso=None)¶ Energy Minimizer (FIRE).
Parameters: - group (
hoomd.group
) – Particle group to apply minimization to. Deprecated in version 2.2:hoomd.md.integrate.mode_minimize_fire()
now accepts integration methods, such ashoomd.md.integrate.nve()
andhoomd.md.integrate.nph()
. The functions operate on user-defined groups. If group is defined here, automaticallyhoomd.md.integrate.nve()
will be used for integration - dt (float) – This is the maximum step size the minimizer is permitted to use. Consider the stability of the system when setting. (in time units)
- Nmin (int) – Number of steps energy change is negative before allowing \(\alpha\) and \(\delta t\) to adapt.
- finc (float) – Factor to increase \(\delta t\) by
- fdec (float) – Factor to decrease \(\delta t\) by
- alpha_start (float) – Initial (and maximum) \(\alpha\)
- falpha (float) – Factor to decrease \(\alpha t\) by
- ftol (float) – force convergence criteria (in units of force over mass)
- wtol (float) – angular momentum convergence criteria (in units of angular momentum)
- Etol (float) – energy convergence criteria (in energy units)
- min_steps (int) – A minimum number of attempts before convergence criteria are considered
- aniso (bool) – Whether to integrate rotational degrees of freedom (bool), default None (autodetect). Added in version 2.2
New in version 2.1.
Changed in version 2.2.
mode_minimize_fire
uses the Fast Inertial Relaxation Engine (FIRE) algorithm to minimize the energy for a group of particles while keeping all other particles fixed. This method is published in Bitzek, et. al., PRL, 2006.At each time step, \(\delta t\), the algorithm uses the NVE Integrator to generate a x, v, and F, and then adjusts v according to
\[\vec{v} = (1-\alpha)\vec{v} + \alpha \hat{F}|\vec{v}|\]where \(\alpha\) and \(\delta t\) are dynamically adaptive quantities. While a current search has been lowering the energy of system for more than \(N_{min}\) steps, \(\alpha\) is decreased by \(\alpha \rightarrow \alpha f_{alpha}\) and \(\delta t\) is increased by \(\delta t \rightarrow max(\delta t \cdot f_{inc}, \delta t_{max})\). If the energy of the system increases (or stays the same), the velocity of the particles is set to 0, \(\alpha \rightarrow \alpha_{start}\) and \(\delta t \rightarrow \delta t \cdot f_{dec}\). Convergence is determined by both the force per particle and the change in energy per particle dropping below ftol and Etol, respectively or
\[\frac{\sum |F|}{N*\sqrt{N_{dof}}} <ftol \;\; and \;\; \Delta \frac{\sum |E|}{N} < Etol\]where N is the number of particles the minimization is acting over (i.e. the group size) Either of the two criterion can be effectively turned off by setting the tolerance to a large number.
If the minimization is acted over a subset of all the particles in the system, the “other” particles will be kept frozen but will still interact with the particles being moved.
Examples:
fire=integrate.mode_minimize_fire(dt=0.05, ftol=1e-2, Etol=1e-7) nve=integrate.nve(group=group.all()) while not(fire.has_converged()): run(100)
Examples:
fire=integrate.mode_minimize_fire(dt=0.05, ftol=1e-2, Etol=1e-7) nph=integrate.nph(group=group.all(),P=0.0,gamma=.5) while not(fire.has_converged()): run(100)
Note
The algorithm requires a base integrator to update the particle position and velocities. Usually this will be either NVE (to minimize energy) or NPH (to minimize energy and relax the box). The quantity minimized is in any case the energy (not the enthalpy or any other quantity).
Note
As a default setting, the algorithm will start with a \(\delta t = \frac{1}{10} \delta t_{max}\) and attempts at least 10 search steps. In practice, it was found that this prevents the simulation from making too aggressive a first step, but also from quitting before having found a good search direction. The minimum number of attempts can be set by the user.
-
get_energy
()¶ Returns the energy after the last iteration of the minimizer
-
has_converged
()¶ Test if the energy minimizer has converged.
Returns: True when the minimizer has converged. Otherwise, return False.
-
reset
()¶ Reset the minimizer to its initial state.
-
restore_state
()¶ Restore the state information from the file used to initialize the simulations
- group (
-
class
hoomd.md.integrate.
mode_standard
(dt, aniso=None)¶ Enables a variety of standard integration methods.
Parameters: - dt (float) – Each time step of the simulation
hoomd.run()
will advance the real time of the system forward by dt (in time units). - aniso (bool) – Whether to integrate rotational degrees of freedom (bool), default None (autodetect).
mode_standard
performs a standard time step integration technique to move the system forward. At each time step, all of the specified forces are evaluated and used in moving the system forward to the next step.By itself,
mode_standard
does nothing. You must specify one or more integration methods to apply to the system. Each integration method can be applied to only a specific group of particles enabling advanced simulation techniques.The following commands can be used to specify the integration methods used by integrate.mode_standard.
There can only be one integration mode active at a time. If there are more than one
integrate.mode_*
commands in a hoomd script, only the most recent before a givenhoomd.run()
will take effect.Examples:
integrate.mode_standard(dt=0.005) integrator_mode = integrate.mode_standard(dt=0.001)
Some integration methods (notable
nvt
,npt
andnph
maintain state between differenthoomd.run()
commands, to allow for restartable simulations. After adding or removing particles, however, a newhoomd.run()
will continue from the old state and the integrator variables will re-equilibrate. To ensure equilibration from a unique reference state (such as all integrator variables set to zero), the method :py:method:reset_methods() can be use to re-initialize the variables.-
reset_methods
()¶ (Re-)initialize the integrator variables in all integration methods
New in version 2.2.
Examples:
run(100) # .. modify the system state, e.g. add particles .. integrator_mode.reset_methods() run(100)
-
restore_state
()¶ Restore the state information from the file used to initialize the simulations
-
set_params
(dt=None, aniso=None)¶ Changes parameters of an existing integration mode.
Parameters: Examples:
integrator_mode.set_params(dt=0.007) integrator_mode.set_params(dt=0.005, aniso=False)
- dt (float) – Each time step of the simulation
-
class
hoomd.md.integrate.
nph
(**params)¶ NPH Integration via MTK barostat-thermostat..
Parameters: nph
performs constant pressure (NPH) simulations using a Martyna-Tobias-Klein barostat, an explicitly reversible and measure-preserving integration scheme. It allows for fully deformable simulation cells and uses the same underlying integrator asnpt
(with nph=True).The available options are identical to those of
npt
, except that kT cannot be specified. For further information, refer to the documentation ofnpt
.Note
A time scale tauP for the relaxation of the barostat is required. This is defined as the relaxation time the barostat would have at an average temperature T_0 = 1, and it is related to the internally used (Andersen) Barostat mass \(W\) via \(W=d N T_0 \tau_P^2\), where \(d\) is the dimensionality and \(N\) the number of particles.
nph
is an integration method and must be used withmode_standard
.Examples:
# Triclinic unit cell nph=integrate.nph(group=all, P=2.0, tauP=1.0, couple="none", all=True) # Cubic unit cell nph = integrate.nph(group=all, P=2.0, tauP=1.0) # Relax the box nph = integrate.nph(group=all, P=0, tauP=1.0, gamma=0.1)
-
disable
()¶ Disables the integration method.
Examples:
method.disable()
Executing the disable command will remove the integration method from the simulation. Any
hoomd.run()
command executed after disabling an integration method will not apply the integration method to the particles during the simulation. A disabled integration method can be re-enabled withenable()
.
-
randomize_velocities
(kT, seed)¶ Assign random velocities and angular momenta to particles in the group, sampling from the Maxwell-Boltzmann distribution. This method considers the dimensionality of the system and particle anisotropy, and removes drift (the center of mass velocity).
New in version 2.3.
Starting in version 2.5, randomize_velocities also chooses random values for the internal integrator variables.
Parameters: Note
Randomization is applied at the start of the next call to
hoomd.run()
.Example:
integrator = md.integrate.nph(group=group.all(), P=2.0, tauP=1.0) integrator.randomize_velocities(kT=1.0, seed=42) run(100)
-
set_params
(kT=None, tau=None, S=None, P=None, tauP=None, rescale_all=None, gamma=None)¶ Changes parameters of an existing integrator.
Parameters: - kT (
hoomd.variant
orfloat
) – New temperature (if set) (in energy units) - tau (float) – New coupling constant (if set) (in time units)
- S (
list
ofhoomd.variant
orfloat
) – New stress components set point (if set) for the barostat (in pressure units). In Voigt notation: [Sxx, Syy, Szz, Syz, Sxz, Sxy] - P (
hoomd.variant
orfloat
) – New isotropic pressure set point (if set) for the barostat (in pressure units). Overrides S if set. - tauP (float) – New barostat coupling constant (if set) (in time units)
- rescale_all (bool) – When True, rescale all particles, not only those in the group
Examples:
integrator.set_params(tau=0.6) integrator.set_params(dt=3e-3, kT=2.0, P=1.0)
- kT (
-
-
class
hoomd.md.integrate.
npt
(group, kT=None, tau=None, S=None, P=None, tauP=None, couple='xyz', x=True, y=True, z=True, xy=False, xz=False, yz=False, all=False, nph=False, rescale_all=None, gamma=None)¶ NPT Integration via MTK barostat-thermostat.
Parameters: - group (
hoomd.group
) – Group of particles on which to apply this method. - kT (
hoomd.variant
orfloat
) – Temperature set point for the thermostat, not needed if nph=True (in energy units). - tau (float) – Coupling constant for the thermostat, not needed if nph=True (in time units).
- S (
list
ofhoomd.variant
orfloat
) – Stress components set point for the barostat (in pressure units). In Voigt notation: [Sxx, Syy, Szz, Syz, Sxz, Sxy] - P (
hoomd.variant
orfloat
) – Isotropic pressure set point for the barostat (in pressure units). Overrides S if set. - tauP (float) – Coupling constant for the barostat (in time units).
- couple (str) – Couplings of diagonal elements of the stress tensor, can be “none”, “xy”, “xz”,”yz”, or “xyz” (default).
- x (bool) – if True, rescale Lx and x component of particle coordinates and velocities
- y (bool) – if True, rescale Ly and y component of particle coordinates and velocities
- z (bool) – if True, rescale Lz and z component of particle coordinates and velocities
- xy (bool) – if True, rescale xy tilt factor and x and y components of particle coordinates and velocities
- xz (bool) – if True, rescale xz tilt factor and x and z components of particle coordinates and velocities
- yz (bool) – if True, rescale yz tilt factor and y and z components of particle coordinates and velocities
- all (bool) – if True, rescale all lengths and tilt factors and components of particle coordinates and velocities
- nph (bool) – if True, integrate without a thermostat, i.e. in the NPH ensemble
- rescale_all (bool) – if True, rescale all particles, not only those in the group
- gamma – (
float
): Dimensionless damping factor for the box degrees of freedom (default: 0)
npt
performs constant pressure, constant temperature simulations, allowing for a fully deformable simulation box.The integration method is based on the rigorous Martyna-Tobias-Klein equations of motion for NPT. For optimal stability, the update equations leave the phase-space measure invariant and are manifestly time-reversible.
By default,
npt
performs integration in a cubic box under hydrostatic pressure by simultaneously rescaling the lengths Lx, Ly and Lz of the simulation box.npt
can also perform more advanced integration modes. The integration mode is specified by a set of couplings and by specifying the box degrees of freedom that are put under barostat control.Couplings define which diagonal elements of the pressure tensor \(P_{\alpha,\beta}\) should be averaged over, so that the corresponding box lengths are rescaled by the same amount.
Valid couplings are:
- none (all box lengths are updated independently)
- xy (Lx and Ly are coupled)
- xz (Lx and Lz are coupled)
- yz (Ly and Lz are coupled)
- xyz (Lx and Ly and Lz are coupled)
The default coupling is xyz, i.e. the ratios between all box lengths stay constant.
Degrees of freedom of the box specify which lengths and tilt factors of the box should be updated, and how particle coordinates and velocities should be rescaled.
Valid keywords for degrees of freedom are:
- x (the box length Lx is updated)
- y (the box length Ly is updated)
- z (the box length Lz is updated)
- xy (the tilt factor xy is updated)
- xz (the tilt factor xz is updated)
- yz (the tilt factor yz is updated)
- all (all elements are updated, equivalent to x, y, z, xy, xz, and yz together)
Any of the six keywords can be combined together. By default, the x, y, and z degrees of freedom are updated.
Note
If any of the diagonal x, y, z degrees of freedom is not being integrated, pressure tensor components along that direction are not considered for the remaining degrees of freedom.
For example:
- Specifying xyz couplings and x, y, and z degrees of freedom amounts to cubic symmetry (default)
- Specifying xy couplings and x, y, and z degrees of freedom amounts to tetragonal symmetry.
- Specifying no couplings and all degrees of freedom amounts to a fully deformable triclinic unit cell
npt
Can also apply a constant stress to the simulation box. To do so, specify the symmetric stress tensor S instead of an isotropic pressure P.Note
npt
assumes that isotropic pressures are positive. Conventions for the stress tensor sometimes assume negative values on the diagonal. You need to set these values negative manually in HOOMD.npt
is an integration method. It must be used withmode_standard
.npt
uses the proper number of degrees of freedom to compute the temperature and pressure of the system in both 2 and 3 dimensional systems, as long as the number of dimensions is set before thenpt
command is specified.For the MTK equations of motion, see:
- G. J. Martyna, D. J. Tobias, M. L. Klein 1994
- M. E. Tuckerman et. al. 2006
- T. Yu et. al. 2010
- Glaser et. al (2013), to be published
Both kT and P can be variant types, allowing for temperature/pressure ramps in simulation runs.
\(\tau\) is related to the Nosé mass \(Q\) by
\[\tau = \sqrt{\frac{Q}{g k_B T_0}}\]where \(g\) is the number of degrees of freedom, and \(k_B T_0\) is the set point (kT above).
A
hoomd.compute.thermo
is automatically specified and associated with group.Examples:
integrate.npt(group=all, kT=1.0, tau=0.5, tauP=1.0, P=2.0) integrator = integrate.npt(group=all, tau=1.0, kT=0.65, tauP = 1.2, P=2.0) # orthorhombic symmetry integrator = integrate.npt(group=all, tau=1.0, kT=0.65, tauP = 1.2, P=2.0, couple="none") # tetragonal symmetry integrator = integrate.npt(group=all, tau=1.0, kT=0.65, tauP = 1.2, P=2.0, couple="xy") # triclinic symmetry integrator = integrate.npt(group=all, tau=1.0, kT=0.65, tauP = 1.2, P=2.0, couple="none", rescale_all=True)
-
disable
()¶ Disables the integration method.
Examples:
method.disable()
Executing the disable command will remove the integration method from the simulation. Any
hoomd.run()
command executed after disabling an integration method will not apply the integration method to the particles during the simulation. A disabled integration method can be re-enabled withenable()
.
-
randomize_velocities
(seed)¶ Assign random velocities and angular momenta to particles in the group, sampling from the Maxwell-Boltzmann distribution. This method considers the dimensionality of the system and particle anisotropy, and removes drift (the center of mass velocity).
New in version 2.3.
Starting in version 2.5, randomize_velocities also chooses random values for the internal integrator variables.
Parameters: seed (int) – Random number seed Note
Randomization is applied at the start of the next call to
hoomd.run()
.Example:
integrator = md.integrate.npt(group=group.all(), kT=1.0, tau=0.5, tauP=1.0, P=2.0) integrator.randomize_velocities(seed=42) run(100)
-
set_params
(kT=None, tau=None, S=None, P=None, tauP=None, rescale_all=None, gamma=None)¶ Changes parameters of an existing integrator.
Parameters: - kT (
hoomd.variant
orfloat
) – New temperature (if set) (in energy units) - tau (float) – New coupling constant (if set) (in time units)
- S (
list
ofhoomd.variant
orfloat
) – New stress components set point (if set) for the barostat (in pressure units). In Voigt notation: [Sxx, Syy, Szz, Syz, Sxz, Sxy] - P (
hoomd.variant
orfloat
) – New isotropic pressure set point (if set) for the barostat (in pressure units). Overrides S if set. - tauP (float) – New barostat coupling constant (if set) (in time units)
- rescale_all (bool) – When True, rescale all particles, not only those in the group
Examples:
integrator.set_params(tau=0.6) integrator.set_params(dt=3e-3, kT=2.0, P=1.0)
- kT (
- group (
-
class
hoomd.md.integrate.
nve
(group, limit=None, zero_force=False)¶ NVE Integration via Velocity-Verlet
Parameters: - group (
hoomd.group
) – Group of particles on which to apply this method. - limit (bool) – (optional) Enforce that no particle moves more than a distance of a limit in a single time step
- zero_force (bool) – When set to true, particles in the a group are integrated forward in time with constant velocity and any net force on them is ignored.
nve
performs constant volume, constant energy simulations using the standard Velocity-Verlet method. For poor initial conditions that include overlapping atoms, a limit can be specified to the movement a particle is allowed to make in one time step. After a few thousand time steps with the limit set, the system should be in a safe state to continue with unconstrained integration.Another use-case for
nve
is to fix the velocity of a certain group of particles. This can be achieved by setting the velocity of those particles in the initial condition and setting the zero_force option to True for that group. A True value for zero_force causes integrate.nve to ignore any net force on each particle and integrate them forward in time with a constant velocity.Note
With an active limit, Newton’s third law is effectively not obeyed and the system can gain linear momentum. Activate the
hoomd.md.update.zero_momentum
updater during the limited nve run to prevent this.nve
is an integration method. It must be used withmode_standard
.A
hoomd.compute.thermo
is automatically specified and associated with group.Examples:
all = group.all() integrate.nve(group=all) integrator = integrate.nve(group=all) typeA = group.type('A') integrate.nve(group=typeA, limit=0.01) integrate.nve(group=typeA, zero_force=True)
-
disable
()¶ Disables the integration method.
Examples:
method.disable()
Executing the disable command will remove the integration method from the simulation. Any
hoomd.run()
command executed after disabling an integration method will not apply the integration method to the particles during the simulation. A disabled integration method can be re-enabled withenable()
.
-
randomize_velocities
(kT, seed)¶ Assign random velocities and angular momenta to particles in the group, sampling from the Maxwell-Boltzmann distribution. This method considers the dimensionality of the system and particle anisotropy, and removes drift (the center of mass velocity).
New in version 2.3.
Parameters: Note
Randomization is applied at the start of the next call to
hoomd.run()
.Example:
integrator = md.integrate.nve(group=group.all()) integrator.randomize_velocities(kT=1.0, seed=42) run(100)
-
set_params
(limit=None, zero_force=None)¶ Changes parameters of an existing integrator.
Parameters: Examples:
integrator.set_params(limit=0.01) integrator.set_params(limit=False)
- group (
-
class
hoomd.md.integrate.
nvt
(group, kT, tau)¶ NVT Integration via the Nosé-Hoover thermostat.
Parameters: - group (
hoomd.group
) – Group of particles on which to apply this method. - kT (
hoomd.variant
orfloat
) – Temperature set point for the Nosé-Hoover thermostat. (in energy units). - tau (float) – Coupling constant for the Nosé-Hoover thermostat. (in time units).
nvt
performs constant volume, constant temperature simulations using the Nosé-Hoover thermostat, using the MTK equations described in Refs. G. J. Martyna, D. J. Tobias, M. L. Klein 1994 and J. Cao, G. J. Martyna 1996.nvt
is an integration method. It must be used in connection withmode_standard
.nvt
uses the proper number of degrees of freedom to compute the temperature of the system in both 2 and 3 dimensional systems, as long as the number of dimensions is set before the integrate.nvt command is specified.\(\tau\) is related to the Nosé mass \(Q\) by
\[\tau = \sqrt{\frac{Q}{g k_B T_0}}\]where \(g\) is the number of degrees of freedom, and \(k_B T_0\) is the set point (kT above).
kT can be a variant type, allowing for temperature ramps in simulation runs.
A
hoomd.compute.thermo
is automatically specified and associated with group.Examples:
all = group.all() integrate.nvt(group=all, kT=1.0, tau=0.5) integrator = integrate.nvt(group=all, tau=1.0, kT=0.65) typeA = group.type('A') integrator = integrate.nvt(group=typeA, tau=1.0, kT=hoomd.variant.linear_interp([(0, 4.0), (1e6, 1.0)]))
-
disable
()¶ Disables the integration method.
Examples:
method.disable()
Executing the disable command will remove the integration method from the simulation. Any
hoomd.run()
command executed after disabling an integration method will not apply the integration method to the particles during the simulation. A disabled integration method can be re-enabled withenable()
.
-
randomize_velocities
(seed)¶ Assign random velocities and angular momenta to particles in the group, sampling from the Maxwell-Boltzmann distribution. This method considers the dimensionality of the system and particle anisotropy, and removes drift (the center of mass velocity).
New in version 2.3.
Starting in version 2.5, randomize_velocities also chooses random values for the internal integrator variables.
Parameters: seed (int) – Random number seed Note
Randomization is applied at the start of the next call to
hoomd.run()
.Example:
integrator = md.integrate.nvt(group=group.all(), kT=1.0, tau=0.5) integrator.randomize_velocities(seed=42) run(100)
- group (