Nose-Hoover Thermostat and Barostat

Overview

Implements Nose-Hoover dynamics using the equations of motion as outlined in Tuckerman et. al, J. Phys. A: Math. Gen 39 (2006) 5629-5651. Allows for dynamics from the NVT or NPT ensembles to be simulated. The total Liouville operator implemented for NPT dynamics is given by (from Tuckerman et. al 2006, p.5641):

\[iL = iL_1 + iL_2 + iL_{\epsilon,1} + iL_{\epsilon,2} + iL_{T-baro} + iL_{T-Part}\]\[iL_1 = \sum\limits_{i=1}^N \bigl[\frac{\mathbf{p}_i}{m_i} + \frac{p_{\epsilon}}{W} \mathbf{r}_i \bigl] \cdot \frac{\partial}{\partial \mathbf{r}_i}\]\[iL_2 = \sum\limits_{i=1}^N \bigl[\mathbf{F}_i - \alpha \frac{p_{\epsilon}}{W}\mathbf{p}_i \bigl] \cdot \frac{\partial}{\partial \mathbf{p}_i}\]\[iL_{\epsilon,1} = \frac{p_{\epsilon}}{W} \frac{\partial}{\partial \epsilon}\]\[iL_{\epsilon,2} = G_{\epsilon} \frac{\partial}{\partial p_{\epsilon}}\]\[\text{where } G_{\epsilon} = \alpha \sum\limits_i \frac{\mathbf{p}_i^2}{m_i} + \sum\limits_{i=1}^N \mathbf{r}_i \cdot \mathbf{F}_i - 3 V \frac{\partial U}{\partial V} - PV\]

Here, \(\mathbf{p}_i\) and \(\mathbf{r}_i\) are the particle momenta and positions, \(\mathbf{F}_i\) are the forces on the particles, \(p_{\epsilon}\) and \(W\) are the barostat momenta and masses, \(\alpha\) is a factor of \(1+\frac{1}{N}\), and \(P\) and \(V\) are the set point pressure and instantaneous volume, respectively.

The corresponding propagator for the NPT ensemble is then given by:

\[\begin{split}\exp(iL \Delta t) = \exp (iL_{T-baro} \frac{\Delta t}{2}) \exp (iL_{T-part} \frac{\Delta t}{2}) \exp (iL_{\epsilon,2} \frac{\Delta t}{2}) \\ \times \exp (iL_2 \frac{\Delta t}{2}) \exp (iL_{\epsilon,1} \Delta t) \exp(iL_1 \Delta t) \exp(iL_2 \frac{\Delta t}{2}) \\ \times \exp(iL_{\epsilon,2} \frac{\Delta t}{2}) \exp(iL_{T-part} \frac{\Delta t}{2}) \exp(iL_{T-baro} \frac{\Delta t}{2})\end{split}\]

The barostat variables and particles are separately thermostatted; in each case, a chain of 3 thermostats is used. Integration is accomplished via a first order Suzuki-Yoshida integration scheme.

Constructor

FixNoseHoover(state,handle,groupHandle)

Arguments

state

The simulation state to which the fix is applied.

handle

The name of the fix. String type.

groupHandle

The group of atoms to which this fix is to be applied. String type.

Python Member Functions

The Nose-Hoover Barostat/Thermostat set points are set via the Python member functions.

Temperature set points may be input through any of the following commands:

setTemperature(temp,timeConstant)
setTemperature(temps,intervals,timeConstant)
setTemperature(tempFunc, timeConstant)

The first method holds the temperature of groupHandle at a set point of temp. The second linearly interpolates between the temperatures given in the list temps at the fractions through the current run given by intervals. The third allows the set point to be determined at each turn via tempFunc. The same scheme is used to set pressure.

Arguments:

temp

The temperature set point for the simulation. Float type.

timeConstant

The time constant associated with the thermostat variables. Float type.

temps

A list of temperature set points to be used throughout the simulation. List of floats.

intervals

A list of fractions through the next run which correspond to the temperatures given in temps. List of floats.

tempFunc

The temperature set point, implemented as a python function.

Likewise, specification of a set point pressure may be accomplished through any of the following commands:

setPressure(mode,pressure,timeConstant)
setPressure(mode,pressures,intervals,timeConstant)
setPressure(mode,pressFunc,timeConstant)

Arguments:

mode

The mode in which cell deformations occur; options are “ISO” or “ANISO”. With mode “ISO”, the internal stress tensor is averaged across the three normal components (or 2, for 2D simulations), and a uniform scale factor for the dimensions emerges. For “ANISO”, the components of the internal stress tensor are not averaged and the individual dimensions are scaled independently.

pressure

The set point pressure for the simulation. Float type.

timeConstant

The time constant associated with the barostat variables. Float type.

pressures

A list of pressure set points to be used through the simulation. List of floats.

intervals

A list of fractions through the next run which correspond to the pressures given in pressures. List of floats.

pressFunc

The pressure set point, implemented as a python function.

For NPT dynamics, both setTemperature and setPressure should be called. They can be called in any order before the simulation is run.

Examples

Example 1: Nose-Hoover Thermostat (NVT Ensemble) - constant set point temperature

# create a simulation state
state = State()

# make an instance of the fix
fixNVT = FixNoseHoover(state, "nvt", "all")

# assign a set point temperature of 300K with time constant 100*state.dt
fixNVT.setTemperature(300.0, 100*state.dt)

# activate the fix
state.activateFix(fixNVT)

Example 2: Nose-Hoover Barostat & Thermostat (NPT Ensemble) - constant set point temperature & pressure

# create a simulation state
state = State()

# make an instance of the fix
fixNPT = FixNoseHoover(state, "npt", "all")

# assign a set point temperature and time constant 100*state.dt
fixNPT.setTemperature(250.0, 100*state.dt)

# assign a set point pressure and time constant 1000*state.dt with isotropic cell deformations
fixNPT.setPressure("ISO", 1.0, 1000*state.dt)

# activate the fix
state.activateFix(fixNPT)

Example 3: Setting temperature via temperature and intervals

# create a simulation state
state = State()

# make an instance of the fix
fixNVT = FixNoseHoover(state, "nvt", "all")

# assign a set point temperature of 300K with time constant 100*state.dt
fixNVT.setTemperature(temps=[100, 500, 400], intervals=[0, 0.2, 1.0], 100*state.dt)

# activate the fix
state.activateFix(fixNVT)

integrator = IntegratorVerlet(state)

#Will sweep from temp=100 to temp=500 between turns 0 and 2000, then 500 and 400 between turns 2000 and 10000
integrator.run(10000)

Example 4: Setting temperature via tempFunc

def randomTemp(turnRunBegan, turnRunEnds, currentTurn):
    return 100 + 50*random() * (currentTurn-turnRunBegan)/(turnRunEnds-turnRunBegan)

# create a simulation state
state = State()

# make an instance of the fix
fixNVT = FixNoseHoover(state, "nvt", "all")

# assign a set point temperature of 300K with time constant 100*state.dt
fixNVT.setTemperature(tempFunc=randomTemp, 100*state.dt)

# activate the fix
state.activateFix(fixNVT)

integrator = IntegratorVerlet(state)

#Will use the value returned by randomTemp each turn as the setpoint
integrator.run(10000)