| 1432 |
|
|
| 1433 |
|
To carry out isobaric-isothermal ensemble calculations {\sc oopse} |
| 1434 |
|
implements the Melchionna modifications to the Nos\'e-Hoover-Andersen |
| 1435 |
< |
equations of motion,\cite{melchionna93} |
| 1435 |
> |
equations of motion.\cite{melchionna93} The equations of motion are the same as NVT with the following exceptions: |
| 1436 |
|
|
| 1437 |
|
\begin{eqnarray} |
| 1438 |
|
\dot{{\bf r}} & = & {\bf v} + \eta \left( {\bf r} - {\bf R}_0 \right), \\ |
| 1439 |
|
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - (\eta + \chi) {\bf v}, \\ |
| 1440 |
– |
\dot{\mathsf{A}} & = & \mathsf{A} \cdot |
| 1441 |
– |
\mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right),\\ |
| 1442 |
– |
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{I}^{-1} |
| 1443 |
– |
\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \cdot \frac{\partial |
| 1444 |
– |
V}{\partial \mathsf{A}} \right) - \chi {\bf j}, \\ |
| 1445 |
– |
\dot{\chi} & = & \frac{1}{\tau_{T}^2} \left( |
| 1446 |
– |
\frac{T}{T_{\mathrm{target}}} - 1 \right) ,\\ |
| 1440 |
|
\dot{\eta} & = & \frac{1}{\tau_{B}^2 f k_B T_{\mathrm{target}}} V \left( P - |
| 1441 |
|
P_{\mathrm{target}} \right), \\ |
| 1442 |
|
\dot{\mathcal{V}} & = & 3 \mathcal{V} \eta . |
| 1485 |
|
file. The units for {\tt tauBarostat} are fs, and the units for the |
| 1486 |
|
{\tt targetPressure} are atmospheres. Like in the NVT integrator, the |
| 1487 |
|
integration of the equations of motion is carried out in a |
| 1488 |
< |
velocity-Verlet style 2 part algorithm: |
| 1488 |
> |
velocity-Verlet style 2 part algorithm with only the following differences: |
| 1489 |
|
|
| 1490 |
|
{\tt moveA:} |
| 1491 |
|
\begin{align*} |
| 1499 |
– |
T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} ,\\ |
| 1500 |
– |
% |
| 1492 |
|
P(t) &\leftarrow \left\{{\bf r}(t)\right\}, \left\{{\bf v}(t)\right\} ,\\ |
| 1493 |
|
% |
| 1494 |
|
{\bf v}\left(t + h / 2\right) &\leftarrow {\bf v}(t) |
| 1495 |
|
+ \frac{h}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t) |
| 1496 |
|
\left(\chi(t) + \eta(t) \right) \right), \\ |
| 1506 |
– |
% |
| 1507 |
– |
{\bf j}\left(t + h / 2 \right) &\leftarrow {\bf j}(t) |
| 1508 |
– |
+ \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t) |
| 1509 |
– |
\chi(t) \right), \\ |
| 1510 |
– |
% |
| 1511 |
– |
\mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left(h * |
| 1512 |
– |
{\bf j}(t + h / 2) \overleftrightarrow{\mathsf{I}}^{-1} |
| 1513 |
– |
\right) ,\\ |
| 1514 |
– |
% |
| 1515 |
– |
\chi\left(t + h / 2 \right) &\leftarrow \chi(t) + |
| 1516 |
– |
\frac{h}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1 |
| 1517 |
– |
\right) ,\\ |
| 1497 |
|
% |
| 1498 |
|
\eta(t + h / 2) &\leftarrow \eta(t) + \frac{h |
| 1499 |
|
\mathcal{V}(t)}{2 N k_B T(t) \tau_B^2} \left( P(t) |
| 1508 |
|
\mathsf{H}(t). |
| 1509 |
|
\end{align*} |
| 1510 |
|
|
| 1511 |
< |
Most of these equations are identical to their counterparts in the NVT |
| 1533 |
< |
integrator, but the propagation of positions to time $t + h$ |
| 1511 |
> |
The propagation of positions to time $t + h$ |
| 1512 |
|
depends on the positions at the same time. {\sc oopse} carries out |
| 1513 |
|
this step iteratively (with a limit of 5 passes through the iterative |
| 1514 |
|
loop). Also, the simulation box $\mathsf{H}$ is scaled uniformly for |
| 1528 |
|
|
| 1529 |
|
{\tt moveB:} |
| 1530 |
|
\begin{align*} |
| 1553 |
– |
T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\}, |
| 1554 |
– |
\left\{{\bf j}(t + h)\right\} ,\\ |
| 1555 |
– |
% |
| 1531 |
|
P(t + h) &\leftarrow \left\{{\bf r}(t + h)\right\}, |
| 1532 |
|
\left\{{\bf v}(t + h)\right\}, \\ |
| 1533 |
|
% |
| 1559 |
– |
\chi\left(t + h \right) &\leftarrow \chi\left(t + h / |
| 1560 |
– |
2 \right) + \frac{h}{2 \tau_T^2} \left( \frac{T(t+h)} |
| 1561 |
– |
{T_{\mathrm{target}}} - 1 \right), \\ |
| 1562 |
– |
% |
| 1534 |
|
\eta(t + h) &\leftarrow \eta(t + h / 2) + |
| 1535 |
|
\frac{h \mathcal{V}(t + h)}{2 N k_B T(t + h) |
| 1536 |
|
\tau_B^2} \left( P(t + h) - P_{\mathrm{target}} \right), \\ |
| 1587 |
|
{\it shape} as well as in the volume of the box. This method utilizes |
| 1588 |
|
the full $3 \times 3$ pressure tensor and introduces a tensor of |
| 1589 |
|
extended variables ($\overleftrightarrow{\eta}$) to control changes to |
| 1590 |
< |
the box shape. The equations of motion for this method are |
| 1590 |
> |
the box shape. The equations of motion for this method differ from those of NPTi as follows: |
| 1591 |
|
\begin{eqnarray} |
| 1592 |
|
\dot{{\bf r}} & = & {\bf v} + \overleftrightarrow{\eta} \cdot \left( {\bf r} - {\bf R}_0 \right), \\ |
| 1593 |
|
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - (\overleftrightarrow{\eta} + |
| 1594 |
|
\chi \cdot \mathsf{1}) {\bf v}, \\ |
| 1624 |
– |
\dot{\mathsf{A}} & = & \mathsf{A} \cdot |
| 1625 |
– |
\mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right) ,\\ |
| 1626 |
– |
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{I}^{-1} |
| 1627 |
– |
\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \cdot \frac{\partial |
| 1628 |
– |
V}{\partial \mathsf{A}} \right) - \chi {\bf j} ,\\ |
| 1629 |
– |
\dot{\chi} & = & \frac{1}{\tau_{T}^2} \left( |
| 1630 |
– |
\frac{T}{T_{\mathrm{target}}} - 1 \right) ,\\ |
| 1595 |
|
\dot{\overleftrightarrow{\eta}} & = & \frac{1}{\tau_{B}^2 f k_B |
| 1596 |
|
T_{\mathrm{target}}} V \left( \overleftrightarrow{\mathsf{P}} - P_{\mathrm{target}}\mathsf{1} \right) ,\\ |
| 1597 |
|
\dot{\mathsf{H}} & = & \overleftrightarrow{\eta} \cdot \mathsf{H} . |
| 1607 |
|
|
| 1608 |
|
{\tt moveA:} |
| 1609 |
|
\begin{align*} |
| 1646 |
– |
T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} ,\\ |
| 1647 |
– |
% |
| 1610 |
|
\overleftrightarrow{\mathsf{P}}(t) &\leftarrow \left\{{\bf r}(t)\right\}, |
| 1611 |
|
\left\{{\bf v}(t)\right\} ,\\ |
| 1612 |
|
% |
| 1615 |
|
\left(\chi(t)\mathsf{1} + \overleftrightarrow{\eta}(t) \right) \cdot |
| 1616 |
|
{\bf v}(t) \right), \\ |
| 1617 |
|
% |
| 1656 |
– |
{\bf j}\left(t + h / 2 \right) &\leftarrow {\bf j}(t) |
| 1657 |
– |
+ \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t) |
| 1658 |
– |
\chi(t) \right), \\ |
| 1659 |
– |
% |
| 1660 |
– |
\mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left(h * |
| 1661 |
– |
{\bf j}(t + h / 2) \overleftrightarrow{\mathsf{I}}^{-1} |
| 1662 |
– |
\right), \\ |
| 1663 |
– |
% |
| 1664 |
– |
\chi\left(t + h / 2 \right) &\leftarrow \chi(t) + |
| 1665 |
– |
\frac{h}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} |
| 1666 |
– |
- 1 \right), \\ |
| 1667 |
– |
% |
| 1618 |
|
\overleftrightarrow{\eta}(t + h / 2) &\leftarrow |
| 1619 |
|
\overleftrightarrow{\eta}(t) + \frac{h \mathcal{V}(t)}{2 N k_B |
| 1620 |
|
T(t) \tau_B^2} \left( \overleftrightarrow{\mathsf{P}}(t) |
| 1636 |
|
|
| 1637 |
|
{\tt moveB:} |
| 1638 |
|
\begin{align*} |
| 1689 |
– |
T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\}, |
| 1690 |
– |
\left\{{\bf j}(t + h)\right\}, \\ |
| 1691 |
– |
% |
| 1639 |
|
\overleftrightarrow{\mathsf{P}}(t + h) &\leftarrow \left\{{\bf r} |
| 1640 |
|
(t + h)\right\}, \left\{{\bf v}(t |
| 1641 |
|
+ h)\right\}, \left\{{\bf f}(t + h)\right\} ,\\ |
| 1642 |
|
% |
| 1696 |
– |
\chi\left(t + h \right) &\leftarrow \chi\left(t + h / |
| 1697 |
– |
2 \right) + \frac{h}{2 \tau_T^2} \left( \frac{T(t+ |
| 1698 |
– |
h)}{T_{\mathrm{target}}} - 1 \right), \\ |
| 1699 |
– |
% |
| 1643 |
|
\overleftrightarrow{\eta}(t + h) &\leftarrow |
| 1644 |
|
\overleftrightarrow{\eta}(t + h / 2) + |
| 1645 |
|
\frac{h \mathcal{V}(t + h)}{2 N k_B T(t + h) |
| 1651 |
|
\frac{{\bf f}(t + h)}{m} - |
| 1652 |
|
(\chi(t + h)\mathsf{1} + \overleftrightarrow{\eta}(t |
| 1653 |
|
+ h)) \right) \cdot {\bf v}(t + h), \\ |
| 1711 |
– |
% |
| 1712 |
– |
{\bf j}\left(t + h \right) &\leftarrow {\bf j}\left(t |
| 1713 |
– |
+ h / 2 \right) + \frac{h}{2} \left( {\bf \tau}^b(t |
| 1714 |
– |
+ h) - {\bf j}(t + h) \chi(t + h) \right) . |
| 1654 |
|
\end{align*} |
| 1655 |
|
|
| 1656 |
|
The iterative schemes for both {\tt moveA} and {\tt moveB} are |
| 1792 |
|
F_{z_{\text{Harmonic}}}(t)=-\frac{\partial U(t)}{\partial z(t)}= |
| 1793 |
|
-k_{\text{Harmonic}}(z(t)-z_{\text{cons}}). |
| 1794 |
|
\end{equation} |
| 1795 |
+ |
Parameters concerning the z-constraint method are summarized in Table~\ref{table:zconParams}. |
| 1796 |
|
|
| 1797 |
+ |
\begin{table} |
| 1798 |
+ |
\caption{The Global Keywords: Z-Constraint Parameters} |
| 1799 |
+ |
\label{table:zconParams} |
| 1800 |
+ |
\begin{center} |
| 1801 |
+ |
% Note when adding or removing columns, the \hsize numbers must add up to the total number |
| 1802 |
+ |
% of columns. |
| 1803 |
+ |
\begin{tabularx}{\linewidth}% |
| 1804 |
+ |
{>{\setlength{\hsize}{1.00\hsize}}X% |
| 1805 |
+ |
>{\setlength{\hsize}{0.4\hsize}}X% |
| 1806 |
+ |
>{\setlength{\hsize}{1.2\hsize}}X% |
| 1807 |
+ |
>{\setlength{\hsize}{1.4\hsize}}X} |
| 1808 |
+ |
|
| 1809 |
+ |
{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline |
| 1810 |
+ |
|
| 1811 |
+ |
{\tt zconsTime} & fs & Sets the frequency at which the {\tt .fz} file is written. & Default sets the frequency to the {\tt runTime} \\ |
| 1812 |
+ |
{\tt nZconstraints} & integer & The number of zconstraint molecules& If using zconstraint method, {\tt nZconstraints} must be set \\ |
| 1813 |
+ |
{\tt zconsForcePolicy} & string& The strategy of subtracting zconstraint force from unconstraint molecules & Possible strategies are BYMASS and BYNUMBER. Default strategy is set to BYMASS\\ |
| 1814 |
+ |
{\tt zconsGap} & \r(A) & Set the distance between two adjacent constraint positions& Used mainly in constraining molecules sequentially \\ |
| 1815 |
+ |
{\tt zconsFixtime} & fs & Sets how long the zconstraint molecule is fixed & {\tt zconsGap} must be set if {\tt zconsGap} is already set.\\ |
| 1816 |
+ |
{\tt zconsUsingSMD} &logical & Flag of using Steered Molecular Dynamics or Harmonic Force to move the molecule & Using harmonic force by default\\ |
| 1817 |
+ |
|
| 1818 |
+ |
\end{tabularx} |
| 1819 |
+ |
\end{center} |
| 1820 |
+ |
\end{table} |
| 1821 |
+ |
|
| 1822 |
+ |
|
| 1823 |
+ |
|
| 1824 |
+ |
\section{\label{sec:minimize}Energy Minimization} |
| 1825 |
+ |
|
| 1826 |
+ |
|
| 1827 |
+ |
As one of the basic procedures of molecular modeling, energy minimization |
| 1828 |
+ |
method is used to identify configurations that are stable points on the energy |
| 1829 |
+ |
surface by adjusting the atomic coordinates. Given a potential energy function |
| 1830 |
+ |
$V$ which depends on a set of coordinates, energy minimization algorithm is |
| 1831 |
+ |
developed to find its minimun value. Different from other packages, the |
| 1832 |
+ |
coordinates in OOPSE not only include cartesian coordinates but also euler |
| 1833 |
+ |
angle if directional atom or rigidbody is involved. Unfortunately, due to the |
| 1834 |
+ |
number of local minima and the cost of computation, in most cases, it is |
| 1835 |
+ |
always impossible to identify the global minimum. OOPSE provides two |
| 1836 |
+ |
frequently used first-derivative algorithms, steepest descents and conjugate |
| 1837 |
+ |
gradient, to find a reasonable local minima. |
| 1838 |
+ |
|
| 1839 |
+ |
Given a coordinate set $x_{k}$ and a search direction $d_{k}$, a line search |
| 1840 |
+ |
algorithm is performed along $d_{k}$ to produce $x_{k+1}=x_{k}+$ $\lambda |
| 1841 |
+ |
_{k}d_{k}$. |
| 1842 |
+ |
|
| 1843 |
+ |
In steepest descent algorithm,% |
| 1844 |
+ |
|
| 1845 |
+ |
\begin{equation} |
| 1846 |
+ |
d_{k}=-\nabla V(x_{k}) |
| 1847 |
+ |
\end{equation} |
| 1848 |
+ |
|
| 1849 |
+ |
|
| 1850 |
+ |
Therefore, the gradient and the direction of next step are always orthogonal |
| 1851 |
+ |
which may causes oscillatory behavior in narrow valleys. To overcome this |
| 1852 |
+ |
problem, the Fletcher-Reeves variant of the conjugate algorithm generates |
| 1853 |
+ |
$d_{k+1}$ from the simple recursion% |
| 1854 |
+ |
|
| 1855 |
+ |
\begin{align} |
| 1856 |
+ |
d_{k+1} & =-\nabla V(x_{k+1})+\gamma_{k}d_{k}\\ |
| 1857 |
+ |
\gamma_{k} & =\frac{\nabla V(x_{k+1})^{T}\nabla V(x_{k+1})}{\nabla |
| 1858 |
+ |
V(x_{k})^{T}\nabla V(x_{k})}% |
| 1859 |
+ |
\end{align} |
| 1860 |
+ |
|
| 1861 |
+ |
|
| 1862 |
+ |
The Polak-Ribiere variant of conjugate gradient defines as% |
| 1863 |
+ |
|
| 1864 |
+ |
\begin{equation} |
| 1865 |
+ |
\gamma_{k}=\frac{[\nabla V(x_{k+1})-\nabla V(x)]^{T}\nabla V(x_{k+1})}{\nabla |
| 1866 |
+ |
V(x_{k})^{T}\nabla V(x_{k})}% |
| 1867 |
+ |
\end{equation} |
| 1868 |
+ |
|
| 1869 |
+ |
|
| 1870 |
+ |
The conjugate gradient method assumes that the conformation is close enough to |
| 1871 |
+ |
a local minimum that the potential energy surface is very nearly quadratic. |
| 1872 |
+ |
When initial structure is far from the minimimum, the steepest descent method |
| 1873 |
+ |
can be superiror to conjugate gradient. Hence, steepest descents may generally |
| 1874 |
+ |
be used for the first 10-100 steps of minimization. Another useful feature of |
| 1875 |
+ |
minimization methods in OOPSE is that a modified SHAKE algorithm can be |
| 1876 |
+ |
applied duing the minimization to constraint the bond length. {\tt bass} parameters concerning the minimizer are given in Table~\ref{table:minimizeParams} |
| 1877 |
+ |
|
| 1878 |
+ |
\begin{table} |
| 1879 |
+ |
\caption{The Global Keywords: Energy Minimizer Parameters} |
| 1880 |
+ |
\label{table:minimizeParams} |
| 1881 |
+ |
\begin{center} |
| 1882 |
+ |
% Note when adding or removing columns, the \hsize numbers must add up to the total number |
| 1883 |
+ |
% of columns. |
| 1884 |
+ |
\begin{tabularx}{\linewidth}% |
| 1885 |
+ |
{>{\setlength{\hsize}{1.00\hsize}}X% |
| 1886 |
+ |
>{\setlength{\hsize}{0.4\hsize}}X% |
| 1887 |
+ |
>{\setlength{\hsize}{1.2\hsize}}X% |
| 1888 |
+ |
>{\setlength{\hsize}{1.4\hsize}}X} |
| 1889 |
+ |
|
| 1890 |
+ |
{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline |
| 1891 |
+ |
|
| 1892 |
+ |
{\tt minimizer} & & & \\ |
| 1893 |
+ |
{\tt minMaxIter} & integer & Sets the maximum iteration in energy minimization & Default value is 200\\ |
| 1894 |
+ |
{\tt minWriteFreq} & interger & Sets the frequency at which the {\tt .dump} and {\tt .stat} files are writtern in energy minimization & \\ |
| 1895 |
+ |
{\tt minStepSize} & double & Set the step size of line search & Default value is 0.01\\ |
| 1896 |
+ |
{\tt minFTol} & double & Sets energy tolerance & Default value is $10^(-8)$\\ |
| 1897 |
+ |
{\tt minGTol} & double & Sets gradient tolerance & Default value is $10^(-8)$\\ |
| 1898 |
+ |
{\tt minLSTol} & double & Sets line search tolerance & Default value is $10^(-8)$\\ |
| 1899 |
+ |
{\tt minLSMaxIter} & integer & Sets the maximum iteration of line searching & Default value is 50\\ |
| 1900 |
+ |
|
| 1901 |
+ |
\end{tabularx} |
| 1902 |
+ |
\end{center} |
| 1903 |
+ |
\end{table} |
| 1904 |
+ |
|
| 1905 |
+ |
|
| 1906 |
|
\section{\label{oopseSec:design}Program Design} |
| 1907 |
|
|
| 1908 |
|
\subsection{\label{sec:architecture} {\sc oopse} Architecture} |
