| 5 |
|
DLM method} |
| 6 |
|
|
| 7 |
|
The default method for integrating the equations of motion in {\sc |
| 8 |
< |
oopse} is carried out using a velocity-Verlet version of the |
| 9 |
< |
symplectic splitting method proposed by Dullweber, Leimkuhler and |
| 10 |
< |
McLachlan (DLM).\cite{Dullweber1997} When there are no directional |
| 11 |
< |
atoms or rigid bodies present in the simulation, this integrator |
| 12 |
< |
becomes the standard velocity-Verlet integrator which is known to |
| 13 |
< |
sample the microcanonical (NVE) ensemble. |
| 8 |
> |
oopse} is a velocity-Verlet version of the symplectic splitting method |
| 9 |
> |
proposed by Dullweber, Leimkuhler and McLachlan |
| 10 |
> |
(DLM).\cite{Dullweber1997} When there are no directional atoms or |
| 11 |
> |
rigid bodies present in the simulation, this integrator becomes the |
| 12 |
> |
standard velocity-Verlet integrator which is known to sample the |
| 13 |
> |
microcanonical (NVE) ensemble.\cite{} |
| 14 |
|
|
| 15 |
|
Previous integration methods for orientational motion have problems |
| 16 |
|
that are avoided in the DLM method. Direct propagation of the Euler |
| 20 |
|
bodies has an orientation near $\theta=0$ or $\theta=\pi$. More |
| 21 |
|
modern quaternion-based integration methods have relatively poor |
| 22 |
|
energy conservation. While quaternions work well for orientational |
| 23 |
< |
motion in other ensembles ensembles, the microcanonical ensemble has a |
| 23 |
> |
motion in other ensembles, the microcanonical ensemble has a |
| 24 |
|
constant energy requirement that is quite sensitive to errors in the |
| 25 |
|
equations of motion. An earlier implementation of {\sc oopse} |
| 26 |
|
utilized quaternions for propagation of rotational motion; however, a |
| 30 |
|
|
| 31 |
|
The key difference in the integration method proposed by Dullweber |
| 32 |
|
\emph{et al.} is that the entire rotation matrix is propagated from |
| 33 |
< |
one time step to the next. In the past, this would not have been a |
| 34 |
< |
feasible option, since that the rotation matrix for a single body has |
| 35 |
< |
nine elements compared to 3 or 4 elements for Euler angles and |
| 36 |
< |
quaternions respectively. System memory has become less costly in |
| 37 |
< |
recent times, and this can be translated into substantial benefits in |
| 38 |
< |
energy conservation. |
| 33 |
> |
one time step to the next. In the past, this would not have been |
| 34 |
> |
feasible, since that the rotation matrix for a single body has nine |
| 35 |
> |
elements compared to 3 or 4 elements for Euler angles and quaternions |
| 36 |
> |
respectively. System memory has become less costly in recent times, |
| 37 |
> |
and this can be translated into substantial benefits in energy |
| 38 |
> |
conservation. |
| 39 |
|
|
| 40 |
|
The basic equations of motion being integrated are derived from the |
| 41 |
|
Hamiltonian for conservative systems containing rigid bodies, |
| 42 |
|
\begin{equation} |
| 43 |
< |
H = \sum_{i} \frac{\vec{v}_i^T M_i \vec{v}_i}{2} + \sum_i |
| 44 |
< |
\frac{\vec{j}_i^T I_i \vec{j_i}}{2} + |
| 45 |
< |
V(\left\{\vec{r}_i\right\}, \left\{\vec{\Omega}_i\right\}) |
| 43 |
> |
H = \sum_{i} \frac{{\bf v}_i^T m_i {\bf v}_i}{2} + \sum_i |
| 44 |
> |
\frac{{\bf j}_i^T \cdot {\bf j}_i}{2 \overleftrightarrow{\mathsf{I}}_i} + |
| 45 |
> |
V\left(\left\{{\bf r}\right\}, \left\{\mathsf{A}\right\}\right) |
| 46 |
|
\end{equation} |
| 47 |
< |
where $\vec{r}_i$, $\vec{v}_i$, $\vec{j}_i$ and $I_i$ are the |
| 48 |
< |
cartesian position vector of the center of mass, its velocity, the |
| 49 |
< |
body-frame angular velocity, and the moment of inertia tensor (also in |
| 50 |
< |
the body frame) of particle $i$, respectively. $V$ is the potential |
| 51 |
< |
energy function and $\vec{\Omega}_i$ is a representation of the |
| 52 |
< |
orientation of particle $i$. The equations of motion for the particle |
| 53 |
< |
centers of mass are derived from Hamilton's equations and are quite |
| 54 |
< |
simple, |
| 47 |
> |
where ${\bf r}_i$ and ${\bf v}_i$ are the cartesian position vector |
| 48 |
> |
and velocity of the center of mass of particle $i$, and ${\bf j}_i$ |
| 49 |
> |
and $\overleftrightarrow{\mathsf{I}}_i$ are the body-fixed angular |
| 50 |
> |
momentum and moment of inertia tensor, respectively. $\mathsf{A}_i$ |
| 51 |
> |
is the $3 \times 3$ rotation matrix describing the instantaneous |
| 52 |
> |
orientation of the particle. $V$ is the potential energy function |
| 53 |
> |
which may depend on both the positions $\left\{{\bf r}\right\}$ and |
| 54 |
> |
orientations $\left\{\mathsf{A}\right\}$ of all particles. The |
| 55 |
> |
equations of motion for the particle centers of mass are derived from |
| 56 |
> |
Hamilton's equations and are quite simple, |
| 57 |
|
\begin{eqnarray} |
| 58 |
< |
\frac{d \vec{r}}{d t} & = & \vec{v}(t) \\ |
| 59 |
< |
\frac{d \vec{v}}{d t} & = & \frac{\vec{f}(t)}{m}, |
| 58 |
> |
\dot{{\bf r}} & = & {\bf v} \\ |
| 59 |
> |
\dot{{\bf v}} & = & \frac{{\bf f}}{m} |
| 60 |
|
\end{eqnarray} |
| 61 |
< |
where $\vec{f}(t)$ is the instantaneous force on the centers of mass, |
| 61 |
> |
where ${\bf f}$ is the instantaneous force on the center of mass |
| 62 |
> |
of the particle, |
| 63 |
|
\begin{equation} |
| 64 |
< |
\vec{f}(t) = \frac{\partial}{\partial |
| 65 |
< |
\vec{r}}V(\left\{\vec{r}_i(t)\right\}, |
| 63 |
< |
\left\{\vec{\Omega}_i(t)\right\}). |
| 64 |
> |
{\bf f} = \frac{\partial}{\partial |
| 65 |
> |
{\bf r}} V(\left\{{\bf r}(t)\right\}, \left\{\mathsf{A}(t)\right\}). |
| 66 |
|
\end{equation} |
| 67 |
|
|
| 68 |
< |
The equations of motion for the orientational degrees of freedom |
| 69 |
< |
require switching between the space-fixed (or laboratory-fixed) |
| 68 |
< |
frame and the body-fixed frame. Forces and torques are most |
| 69 |
< |
conveniently calculated in the space-fixed frame, while the rotational |
| 70 |
< |
equations of motion are simplest in the body-fixed frame, |
| 68 |
> |
|
| 69 |
> |
The equations of motion for the orientational degrees of freedom are |
| 70 |
|
\begin{eqnarray} |
| 71 |
< |
\frac{d j_{x}}{d t} & = & \frac{\tau_x(t)}{I_{xx}} + |
| 72 |
< |
\left(\frac{I_{yy} - I_{zz}}{I_{xx}} \right) j_y j_z \\ |
| 73 |
< |
\frac{d j_{y}}{d t} & = & \frac{\tau_y(t)}{I_{yy}} + |
| 74 |
< |
\left(\frac{I_{zz} - I_{xx}}{I_{yy}} \right) j_z j_x \\ |
| 75 |
< |
\frac{d j_{z}}{d t} & = & \frac{\tau_z(t)}{I_{zz}} + |
| 77 |
< |
\left(\frac{I_{xx} - I_{yy}}{I_{zz}} \right) j_x j_y |
| 71 |
> |
\dot{\mathsf{A}} & = & \mathsf{A} |
| 72 |
> |
\mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right) \\ |
| 73 |
> |
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{I}^{-1} |
| 74 |
> |
\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \frac{\partial |
| 75 |
> |
V}{\partial \mathsf{A}} \right) |
| 76 |
|
\end{eqnarray} |
| 77 |
+ |
In these equations of motion, the $\mbox{skew}$ matrix of a vector |
| 78 |
+ |
${\bf v} = \left( v_1, v_2, v_3 \right)$ is defined: |
| 79 |
+ |
\begin{equation} |
| 80 |
+ |
\mbox{skew}\left( {\bf v} \right) := \left( |
| 81 |
+ |
\begin{array}{ccc} |
| 82 |
+ |
0 & v_3 & - v_2 \\ |
| 83 |
+ |
-v_3 & 0 & v_1 \\ |
| 84 |
+ |
v_2 & -v_1 & 0 |
| 85 |
+ |
\end{array} |
| 86 |
+ |
\right) |
| 87 |
+ |
\end{equation} |
| 88 |
+ |
The $\mbox{rot}$ notation refers to the mapping of the $3 \times 3$ |
| 89 |
+ |
rotation matrix to a vector of orientations by first computing the |
| 90 |
+ |
skew-symmetric part $\left(\mathsf{A} - \mathsf{A}^{T}\right)$ and |
| 91 |
+ |
then associating this with a length 3 vector by inverting the |
| 92 |
+ |
$\mbox{skew}$ function above: |
| 93 |
+ |
\begin{equation} |
| 94 |
+ |
\mbox{rot}\left(\mathsf{A}\right) := \mbox{ skew}^{-1}\left(\mathsf{A} |
| 95 |
+ |
- \mathsf{A}^{T} \right) |
| 96 |
+ |
\end{equation} |
| 97 |
|
|
| 80 |
– |
The DLM method allows for Verlet style integration of both linear and |
| 81 |
– |
angular motion of rigid bodies. |
| 98 |
|
|
| 99 |
|
|
| 84 |
– |
|
| 85 |
– |
|
| 86 |
– |
|
| 100 |
|
In the integration method, the |
| 101 |
|
orientational propagation involves a sequence of matrix evaluations to |
| 102 |
|
update the rotation matrix.\cite{Dullweber1997} These matrix rotations |
| 144 |
|
|
| 145 |
|
\subsection{\label{sec:extended}Extended Systems for other Ensembles} |
| 146 |
|
|
| 147 |
+ |
{\sc oopse} implements a number of extended system integrators for |
| 148 |
+ |
sampling from other ensembles relevant to chemical physics. The |
| 149 |
+ |
integrator can selected with the {\tt ensemble} keyword in the |
| 150 |
+ |
{\tt .bass} file: |
| 151 |
|
|
| 152 |
+ |
\begin{center} |
| 153 |
+ |
\begin{tabular}{lll} |
| 154 |
+ |
{\bf Integrator} & {\bf Ensemble} & {\bf {\tt .bass} line} \\ |
| 155 |
+ |
NVE & microcanonical & {\tt ensemble = ``NVE''; } \\ |
| 156 |
+ |
NVT & canonical & {\tt ensemble = ``NVT''; } \\ |
| 157 |
+ |
NPTi & isobaric-isothermal (with isotropic volume changes) & {\tt |
| 158 |
+ |
ensemble = ``NPTi'';} \\ |
| 159 |
+ |
NPTf & isobaric-isothermal (with changes to box shape) & {\tt |
| 160 |
+ |
ensemble = ``NPTf'';} \\ |
| 161 |
+ |
NPTxyz & approximate isobaric-isothermal & {\tt ensemble = |
| 162 |
+ |
``NPTxyz'';} \\ |
| 163 |
+ |
& (with separate barostats on each box dimension) & |
| 164 |
+ |
\end{tabular} |
| 165 |
+ |
\end{center} |
| 166 |
|
|
| 167 |
+ |
The relatively well-known Nos\'e-Hoover thermostat is implemented in |
| 168 |
+ |
{\sc oopse}'s NVT integrator. This method couples an extra degree of |
| 169 |
+ |
freedom (the thermostat) to the kinetic energy of the system, and has |
| 170 |
+ |
been shown to sample the canonical distribution in the system degrees |
| 171 |
+ |
of freedom while conserving a quantity that is, to within a constant, |
| 172 |
+ |
the Helmholtz free energy. |
| 173 |
|
|
| 174 |
< |
{\sc oopse} implements a |
| 174 |
> |
NPT algorithms attempt to maintain constant pressure in the system by |
| 175 |
> |
coupling the volume of the system to a barostat. {\sc oopse} contains |
| 176 |
> |
three different constant pressure algorithms. The first two, NPTi and |
| 177 |
> |
NPTf have been shown to conserve a quantity that is, to within a |
| 178 |
> |
constant, the Gibbs free energy. The Melchionna modification to the |
| 179 |
> |
Hoover barostat is implemented in both NPTi and NPTf. NPTi allows |
| 180 |
> |
only isotropic changes in the simulation box, while box {\it shape} |
| 181 |
> |
variations are allowed in NPTf. The NPTxyz integrator has {\it not} |
| 182 |
> |
been shown to sample from the isobaric-isothermal ensemble. It is |
| 183 |
> |
useful, however, in that it maintains orthogonality for the axes of |
| 184 |
> |
the simulation box while attempting to equalize pressure along the |
| 185 |
> |
three perpendicular directions in the box. |
| 186 |
|
|
| 187 |
+ |
Each of the extended system integrators requires additional keywords |
| 188 |
+ |
to set target values for the thermodynamic state variables that are |
| 189 |
+ |
being held constant. Keywords are also required to set the |
| 190 |
+ |
characteristic decay times for the dynamics of the extended |
| 191 |
+ |
variables. |
| 192 |
|
|
| 193 |
< |
\subsubsection{\label{sec:noseHooverThermo}Nose-Hoover Thermostatting} |
| 193 |
> |
\begin{tabular}{llll} |
| 194 |
> |
{\bf variable} & {\bf {\tt .bass} keyword} & {\bf units} & {\bf |
| 195 |
> |
default value} \\ |
| 196 |
> |
$T_{\mathrm{target}}$ & {\tt targetTemperature = 300;} & K & none \\ |
| 197 |
> |
$P_{\mathrm{target}}$ & {\tt targetPressure = 1;} & atm & none \\ |
| 198 |
> |
$\tau_T$ & {\tt tauThermostat = 1e3;} & fs & none \\ |
| 199 |
> |
$\tau_B$ & {\tt tauBarostat = 5e3;} & fs & none \\ |
| 200 |
> |
& {\tt resetTime = 200;} & fs & none \\ |
| 201 |
> |
& {\tt useInitialExtendedSystemState = ``true'';} & logical & |
| 202 |
> |
false |
| 203 |
> |
\end{tabular} |
| 204 |
|
|
| 205 |
< |
To mimic the effects of being in a constant temperature ({\sc nvt}) |
| 206 |
< |
ensemble, {\sc oopse} uses the Nose-Hoover extended system |
| 207 |
< |
approach.\cite{Hoover85} In this method, the equations of motion for |
| 208 |
< |
the particle positions and velocities are |
| 205 |
> |
Two additional keywords can be used to either clear the extended |
| 206 |
> |
system variables periodically ({\tt resetTime}), or to maintain the |
| 207 |
> |
state of the extended system variables between simulations ({\tt |
| 208 |
> |
useInitialExtendedSystemState}). More details on these variables |
| 209 |
> |
and their use in the integrators follows below. |
| 210 |
> |
|
| 211 |
> |
\subsubsection{\label{sec:noseHooverThermo}Nos\'{e}-Hoover Thermostatting} |
| 212 |
> |
|
| 213 |
> |
The Nos\'e-Hoover equations of motion are given by\cite{Hoover85} |
| 214 |
|
\begin{eqnarray} |
| 215 |
|
\dot{{\bf r}} & = & {\bf v} \\ |
| 216 |
< |
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - \chi {\bf v} |
| 216 |
> |
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - \chi {\bf v} \\ |
| 217 |
> |
\dot{\mathsf{A}} & = & \\ |
| 218 |
> |
\dot{{\bf j}} & = & - \chi {\bf j} |
| 219 |
|
\label{eq:nosehoovereom} |
| 220 |
|
\end{eqnarray} |
| 221 |
|
|
| 222 |
|
$\chi$ is an ``extra'' variable included in the extended system, and |
| 223 |
|
it is propagated using the first order equation of motion |
| 224 |
|
\begin{equation} |
| 225 |
< |
\dot{\chi} = \frac{1}{\tau_{T}} \left( \frac{T}{T_{target}} - 1 \right) |
| 225 |
> |
\dot{\chi} = \frac{1}{\tau_{T}^2} \left( \frac{T}{T_{\mathrm{target}}} - 1 \right). |
| 226 |
|
\label{eq:nosehooverext} |
| 227 |
|
\end{equation} |
| 158 |
– |
where $T_{target}$ is the target temperature for the simulation, and |
| 159 |
– |
$\tau_{T}$ is a time constant for the thermostat. |
| 228 |
|
|
| 229 |
< |
To select the Nose-Hoover {\sc nvt} ensemble, the {\tt ensemble = NVT;} |
| 230 |
< |
command would be used in the simulation's {\sc bass} file. There is |
| 231 |
< |
some subtlety in choosing values for $\tau_{T}$, and it is usually set |
| 232 |
< |
to values of a few ps. Within a {\sc bass} file, $\tau_{T}$ could be |
| 233 |
< |
set to 1 ps using the {\tt tauThermostat = 1000; } command. |
| 229 |
> |
The instantaneous temperature $T$ is proportional to the total kinetic |
| 230 |
> |
energy (both translational and orientational) and is given by |
| 231 |
> |
\begin{equation} |
| 232 |
> |
T = \frac{2 K}{f k_B} |
| 233 |
> |
\end{equation} |
| 234 |
> |
Here, $f$ is the total number of degrees of freedom in the system, |
| 235 |
> |
\begin{equation} |
| 236 |
> |
f = 3 N + 3 N_{\mathrm{orient}} - N_{\mathrm{constraints}} |
| 237 |
> |
\end{equation} |
| 238 |
> |
and $K$ is the total kinetic energy, |
| 239 |
> |
\begin{equation} |
| 240 |
> |
K = \sum_{i=1}^{N} \frac{1}{2} m_i {\bf v}_i^T \cdot {\bf v}_i + |
| 241 |
> |
\sum_{i=1}^{N_{\mathrm{orient}}} \sum_{\alpha=x,y,z} \frac{{\bf |
| 242 |
> |
j}_{i\alpha}^T \cdot {\bf j}_{i\alpha}}{2 |
| 243 |
> |
\overleftrightarrow{\mathsf{I}}_{i,\alpha \alpha}} |
| 244 |
> |
\end{equation} |
| 245 |
|
|
| 246 |
+ |
In eq.(\ref{eq:nosehooverext}), $\tau_T$ is the time constant for |
| 247 |
+ |
relaxation of the temperature to the target value. To set values for |
| 248 |
+ |
$\tau_T$ or $T_{\mathrm{target}}$ in a simulation, one would use the |
| 249 |
+ |
{\tt tauThermostat} and {\tt targetTemperature} keywords in the {\tt |
| 250 |
+ |
.bass} file. The units for {\tt tauThermostat} are fs, and the units |
| 251 |
+ |
for the {\tt targetTemperature} are degrees K. The integration of |
| 252 |
+ |
the equations of motion is carried out in a velocity-Verlet style 2 |
| 253 |
+ |
part algorithm: |
| 254 |
|
|
| 255 |
+ |
{\tt moveA:} |
| 256 |
+ |
\begin{eqnarray} |
| 257 |
+ |
T(t) & \leftarrow & \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\ |
| 258 |
+ |
{\bf v}\left(t + \delta t / 2\right) & \leftarrow & {\bf |
| 259 |
+ |
v}(t) + \frac{\delta t}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t) |
| 260 |
+ |
\chi(t)\right) \\ |
| 261 |
+ |
{\bf r}(t + \delta t) & \leftarrow & {\bf r}(t) + \delta t {\bf |
| 262 |
+ |
v}\left(t + \delta t / 2 \right) \\ |
| 263 |
+ |
{\bf j}\left(t + \delta t / 2 \right) & \leftarrow & {\bf |
| 264 |
+ |
j}(t) + \frac{\delta t}{2} \left( {\bf \tau}^b(t) - {\bf j}(t) |
| 265 |
+ |
\chi(t) \right) \\ |
| 266 |
+ |
\mathsf{A}(t + \delta t) & \leftarrow & \mathrm{rot}({\bf j}(t + |
| 267 |
+ |
\delta t / 2) \overleftrightarrow{\mathsf{I}}^{b}, |
| 268 |
+ |
\delta t) \\ |
| 269 |
+ |
\chi\left(t + \delta t / 2 \right) & \leftarrow & \chi(t) + |
| 270 |
+ |
\frac{\delta t}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1 |
| 271 |
+ |
\right) |
| 272 |
+ |
\end{eqnarray} |
| 273 |
+ |
|
| 274 |
+ |
Here $\mathrm{rot}( {\bf j} / {\bf I} )$ is the same symplectic Trotter |
| 275 |
+ |
factorization of the three rotation operations that was discussed in |
| 276 |
+ |
the section on the DLM integrator. Note that this operation modifies |
| 277 |
+ |
both the rotation matrix $\mathsf{A}$ and the angular momentum ${\bf |
| 278 |
+ |
j}$. {\tt moveA} propagates velocities by a half time step, and |
| 279 |
+ |
positional degrees of freedom by a full time step. The new positions |
| 280 |
+ |
(and orientations) are then used to calculate a new set of forces and |
| 281 |
+ |
torques. |
| 282 |
+ |
|
| 283 |
+ |
{\tt doForces:} |
| 284 |
+ |
\begin{eqnarray} |
| 285 |
+ |
{\bf f}(t + \delta t) & \leftarrow & - \frac{\partial V}{\partial {\bf |
| 286 |
+ |
r}(t + \delta t)} \\ |
| 287 |
+ |
{\bf \tau}^{s}(t + \delta t) & \leftarrow & {\bf u}(t + \delta t) |
| 288 |
+ |
\times \frac{\partial V}{\partial {\bf u}} \\ |
| 289 |
+ |
{\bf \tau}^{b}(t + \delta t) & \leftarrow & \mathsf{A}(t + \delta t) |
| 290 |
+ |
\cdot {\bf \tau}^s(t + \delta t) |
| 291 |
+ |
\end{eqnarray} |
| 292 |
+ |
|
| 293 |
+ |
Here ${\bf u}$ is a unit vector pointing along the principal axis of |
| 294 |
+ |
the rigid body being propagated, ${\bf \tau}^s$ is the torque in the |
| 295 |
+ |
space-fixed (laboratory) frame, and ${\bf \tau}^b$ is the torque in |
| 296 |
+ |
the body-fixed frame. ${\bf u}$ is automatically calculated when the |
| 297 |
+ |
rotation matrix $\mathsf{A}$ is calculated in {\tt moveA}. |
| 298 |
+ |
|
| 299 |
+ |
Once the forces and torques have been obtained at the new time step, |
| 300 |
+ |
the velocities can be advanced to the same time value. |
| 301 |
+ |
|
| 302 |
+ |
{\tt moveB:} |
| 303 |
+ |
\begin{eqnarray} |
| 304 |
+ |
T(t + \delta t) & \leftarrow & \left\{{\bf v}(t + \delta t)\right\}, |
| 305 |
+ |
\left\{{\bf j}(t + \delta t)\right\} \\ |
| 306 |
+ |
\chi\left(t + \delta t \right) & \leftarrow & \chi\left(t + \delta t / |
| 307 |
+ |
2 \right) + \frac{\delta t}{2 \tau_T^2} \left( \frac{T(t+\delta |
| 308 |
+ |
t)}{T_{\mathrm{target}}} - 1 \right) \\ |
| 309 |
+ |
{\bf v}\left(t + \delta t \right) & \leftarrow & {\bf |
| 310 |
+ |
v}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( |
| 311 |
+ |
\frac{{\bf f}(t + \delta t)}{m} - {\bf v}(t + \delta t) |
| 312 |
+ |
\chi(t \delta t)\right) \\ |
| 313 |
+ |
{\bf j}\left(t + \delta t \right) & \leftarrow & {\bf |
| 314 |
+ |
j}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( {\bf |
| 315 |
+ |
\tau}^b(t + \delta t) - {\bf j}(t + \delta t) |
| 316 |
+ |
\chi(t + \delta t) \right) |
| 317 |
+ |
\end{eqnarray} |
| 318 |
+ |
|
| 319 |
+ |
Since ${\bf v}(t + \delta t)$ and ${\bf j}(t + \delta t)$ are required |
| 320 |
+ |
to caclculate $T(t + \delta t)$ as well as $\chi(t + \delta t)$, the |
| 321 |
+ |
indirectly depend on their own values at time $t + \delta t$. {\tt |
| 322 |
+ |
moveB} is therefore done in an iterative fashion until $\chi(t + |
| 323 |
+ |
\delta t)$ becomes self-consistent. The relative tolerance for the |
| 324 |
+ |
self-consistency check defaults to a value of $\mbox{10}^{-6}$, but |
| 325 |
+ |
{\sc oopse} will terminate the iteration after 4 loops even if the |
| 326 |
+ |
consistency check has not been satisfied. |
| 327 |
+ |
|
| 328 |
+ |
The Nos\'e-Hoover algorithm is known to conserve a Hamiltonian for the |
| 329 |
+ |
extended system that is, to within a constant, identical to the |
| 330 |
+ |
Helmholtz free energy, |
| 331 |
+ |
\begin{equation} |
| 332 |
+ |
H_{\mathrm{NVT}} = V + K + f k_B T_{\mathrm{target}} \left( |
| 333 |
+ |
\frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t\prime) dt\prime |
| 334 |
+ |
\right) |
| 335 |
+ |
\end{equation} |
| 336 |
+ |
Poor choices of $\delta t$ or $\tau_T$ can result in non-conservation |
| 337 |
+ |
of $H_{\mathrm{NVT}}$, so the conserved quantity is maintained in the |
| 338 |
+ |
last column of the {\tt .stat} file to allow checks on the quality of |
| 339 |
+ |
the integration. |
| 340 |
+ |
|
| 341 |
+ |
Bond constraints are applied at the end of both the {\tt moveA} and |
| 342 |
+ |
{\tt moveB} portions of the algorithm. Details on the constraint |
| 343 |
+ |
algorithms are given in section \ref{sec:rattle}. |
| 344 |
+ |
|
| 345 |
+ |
\subsubsection{\label{sec:NPTi}Constant-pressure integration (isotropic box)} |
| 346 |
+ |
|
| 347 |
+ |
|
| 348 |
|
\subsection{\label{Sec:zcons}Z-Constraint Method} |
| 349 |
|
|
| 350 |
|
Based on fluctuatin-dissipation theorem,\bigskip\ force auto-correlation |
