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} |