About OOPSE OOPSE News Features Download CVS Mailing Lists Documentation Authors Funding Results	Constant-pressure integration with isotropic box deformations (NPTi) To carry out isobaric-isothermal ensemble calculations, OOPSE implements the Melchionna modifications to the Nosé-Hoover-Andersen equations of motion.[46] The equations of motion are the same as NVT with the following exceptions: $$\dot{{\bf r}} = {\bf v} + \eta \left( {\bf r} - {\bf R}_0 \right),$$ (4.26) $$\dot{{\bf v}} = \frac{{\bf f}}{m} - (\eta + \chi) {\bf v},$$ (4.27) $$\dot{\eta} = \frac{1}{\tau_{B}^2 f k_B T_{\mathrm{target}}} V \left( P - P_{\mathrm{target}} \right),$$ (4.28) $$\dot{\mathcal{V}} = 3 \mathcal{V} \eta .$$ (4.29) $\chi$ and $\eta$ are the ``extra'' degrees of freedom in the extended system. $\chi$ is a thermostat, and it has the same function as it does in the Nosé-Hoover NVT integrator. $\eta$ is a barostat which controls changes to the volume of the simulation box. ${\bf R}_0$ is the location of the center of mass for the entire system, and $\mathcal{V}$ is the volume of the simulation box. At any time, the volume can be calculated from the determinant of the matrix which describes the box shape: $$\mathcal{V} = \det(\mathsf{H}).$$ (4.30) The NPTi integrator requires an instantaneous pressure. This quantity is calculated via the pressure tensor, $$\overleftrightarrow{\mathsf{P}}(t) = \frac{1}{\mathcal{V}(t)} \le... ...{\bf v}_i(t) \otimes {\bf v}_i(t) \right) + \overleftrightarrow{\mathsf{W}}(t).$$ (4.31) The kinetic contribution to the pressure tensor utilizes the outer product of the velocities, denoted by the $\otimes$ symbol. The stress tensor is calculated from another outer product of the inter-atomic separation vectors ( ${\bf r}_{ij} = {\bf r}_j - {\bf r}_i$ ) with the forces between the same two atoms, $$\overleftrightarrow{\mathsf{W}}(t) = \sum_{i} \sum_{j>i} {\bf r}_{ij}(t) \otimes {\bf f}_{ij}(t).$$ (4.32) In systems containing cutoff groups, the stress tensor is computed between the centers-of-mass of the cutoff groups: $$\overleftrightarrow{\mathsf{W}}(t) = \sum_{a} \sum_{b} {\bf r}_{ab}(t) \otimes {\bf f}_{ab}(t).$$ (4.33) where ${\bf r}_{ab}$ is the distance between the centers of mass, and $${\bf f}_{ab} = s(r_{ab}) \sum_{i \in a} \sum_{j \in b} {\bf f}_{i... ...r}_{ab}}{\vert r_{ab}\vert} \sum_{i \in a} \sum_{j \in b} V_{ij}({\bf r}_{ij}).$$ (4.34) The instantaneous pressure is then simply obtained from the trace of the pressure tensor, $$P(t) = \frac{1}{3} \mathrm{Tr} \left( \overleftrightarrow{\mathsf{P}}(t) \right).$$ (4.35) In eq.(4.29), $\tau_B$ is the time constant for relaxation of the pressure to the target value. To set values for $\tau_B$ or $P_{\mathrm{target}}$ in a simulation, one would use the tauBarostat and targetPressure keywords in the meta-data file. The units for tauBarostat are fs, and the units for the targetPressure are atmospheres. Like in the NVT integrator, the integration of the equations of motion is carried out in a velocity-Verlet style two part algorithm with only the following differences: moveA: $$P(t) \leftarrow \left\{{\bf r}(t)\right\}, \left\{{\bf v}(t)\right\} ,$$ $${\bf v}\left(t + h / 2\right) \leftarrow {\bf v}(t) + \frac{h}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t) \left(\chi(t) + \eta(t) \right) \right),$$ $$\eta(t + h / 2) \leftarrow \eta(t) + \frac{h \mathcal{V}(t)}{2 N k_B T(t) \tau_B^2} \left( P(t) - P_{\mathrm{target}} \right),$$ $${\bf r}(t + h) \leftarrow {\bf r}(t) + h \left\{ {\bf v}\left(t + h / 2 \right) + \eta(t + h / 2)\left[ {\bf r}(t + h) - {\bf R}_0 \right] \right\} ,$$ $$\mathsf{H}(t + h) \leftarrow e^{-h \eta(t + h / 2)} \mathsf{H}(t).$$ The propagation of positions to time $t + h$ depends on the positions at the same time. OOPSE carries out this step iteratively (with a limit of 5 passes through the iterative loop). Also, the simulation box $\mathsf{H}$ is scaled uniformly for one full time step by an exponential factor that depends on the value of $\eta$ at time $t + h / 2$ . Reshaping the box uniformly also scales the volume of the box by $$\mathcal{V}(t + h) \leftarrow e^{ - 3 h \eta(t + h /2)} \times \mathcal{V}(t).$$ (4.36) The doForces step for the NPTi integrator is exactly the same as in both the DLM and NVT integrators. Once the forces and torques have been obtained at the new time step, the velocities can be advanced to the same time value. moveB: $$P(t + h) \leftarrow \left\{{\bf r}(t + h)\right\}, \left\{{\bf v}(t + h)\right\},$$ $$\eta(t + h) \leftarrow \eta(t + h / 2) + \frac{h \mathcal{V}(t + h)}{2 N k_B T(t + h) \tau_B^2} \left( P(t + h) - P_{\mathrm{target}} \right),$$ $${\bf v}\left(t + h \right) \leftarrow {\bf v}\left(t + h / 2 \right) + \frac{h}{2} \left( \frac{{\bf f}(t + h)}{m} - {\bf v}(t + h) (\chi(t + h) + \eta(t + h)) \right) ,$$ $${\bf j}\left(t + h \right) \leftarrow {\bf j}\left(t + h / 2 \right) + \frac{h}{2} \left( {\bf\tau}^b(t + h) - {\bf j}(t + h) \chi(t + h) \right) .$$ Once again, since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required to calculate $T(t + h)$ , $P(t + h)$ , $\chi(t + h)$ , and $\eta(t + h)$ , they indirectly depend on their own values at time $t + h$ . moveB is therefore done in an iterative fashion until $\chi(t + h)$ and $\eta(t + h)$ become self-consistent. The relative tolerance for the self-consistency check defaults to a value of 10$^{-6}$ , but OOPSE will terminate the iteration after 4 loops even if the consistency check has not been satisfied. The Melchionna modification of the Nosé-Hoover-Andersen algorithm is known to conserve a Hamiltonian for the extended system that is, to within a constant, identical to the Gibbs free energy, $$H_{\mathrm{NPTi}} = V + K + f k_B T_{\mathrm{target}} \left( \fra... ..._{0}^{t} \chi(t^\prime) dt^\prime \right) + P_{\mathrm{target}} \mathcal{V}(t).$$ (4.37) Poor choices of $\delta t$ , $\tau_T$ , or $\tau_B$ can result in non-conservation of $H_{\mathrm{NPTi}}$ , so the conserved quantity is maintained in the last column of the .stat file to allow checks on the quality of the integration. It is also known that this algorithm samples the equilibrium distribution for the enthalpy (including contributions for the thermostat and barostat), $$H_{\mathrm{NPTi}} = V + K + \frac{f k_B T_{\mathrm{target}}}{2} \... ...\chi^2 \tau_T^2 + \eta^2 \tau_B^2 \right) + P_{\mathrm{target}} \mathcal{V}(t).$$ (4.38) Bond constraints are applied at the end of both the moveA and moveB portions of the algorithm. Details on the constraint algorithms are given in section 4.7.1.
	Copyright © 2006 - University of Notre Dame
	Updated on January 16, 2006