| 279 |
|
\epsilon_{ij} = \sqrt{\epsilon_{ii} \epsilon_{jj}} |
| 280 |
|
\label{eq:epsilonMix} |
| 281 |
|
\end{equation} |
| 282 |
– |
|
| 283 |
– |
|
| 282 |
|
|
| 283 |
|
\subsection{\label{oopseSec:DUFF}Dipolar Unified-Atom Force Field} |
| 284 |
|
|
| 711 |
|
|
| 712 |
|
\subsection{{\sc bass} and Model Files} |
| 713 |
|
|
| 714 |
< |
Every {\sc oopse} simulation begins with a {\sc bass} file. {\sc bass} |
| 715 |
< |
(\underline{B}izarre \underline{A}tom \underline{S}imulation |
| 716 |
< |
\underline{S}yntax) is a script syntax that is parsed by {\sc oopse} at |
| 717 |
< |
runtime. The {\sc bass} file allows for the user to completely describe the |
| 718 |
< |
system they are to simulate, as well as tailor {\sc oopse}'s behavior during |
| 719 |
< |
the simulation. {\sc bass} files are denoted with the extension |
| 720 |
< |
\texttt{.bass}, an example file is shown in |
| 721 |
< |
Fig.~\ref{fig:bassExample}. |
| 714 |
> |
Every {\sc oopse} simulation begins with a Bizarre Atom Simulation |
| 715 |
> |
Syntax ({\sc bass}) file. {\sc bass} is a script syntax that is parsed |
| 716 |
> |
by {\sc oopse} at runtime. The {\sc bass} file allows for the user to |
| 717 |
> |
completely describe the system they wish to simulate, as well as tailor |
| 718 |
> |
{\sc oopse}'s behavior during the simulation. {\sc bass} files are |
| 719 |
> |
denoted with the extension |
| 720 |
> |
\texttt{.bass}, an example file is shown in |
| 721 |
> |
Scheme~\ref{sch:bassExample}. |
| 722 |
|
|
| 723 |
< |
\begin{figure} |
| 726 |
< |
\centering |
| 727 |
< |
\framebox[\linewidth]{\rule{0cm}{0.75\linewidth}I'm a {\sc bass} file!} |
| 728 |
< |
\caption{Here is an example \texttt{.bass} file} |
| 729 |
< |
\label{fig:bassExample} |
| 730 |
< |
\end{figure} |
| 723 |
> |
\begin{lstlisting}[float,caption={[An example of a complete {\sc bass} file] An example showing a complete {\sc bass} file.},label={sch:bassExample}] |
| 724 |
|
|
| 725 |
+ |
molecule{ |
| 726 |
+ |
name = "Ar"; |
| 727 |
+ |
nAtoms = 1; |
| 728 |
+ |
atom[0]{ |
| 729 |
+ |
type="Ar"; |
| 730 |
+ |
position( 0.0, 0.0, 0.0 ); |
| 731 |
+ |
} |
| 732 |
+ |
} |
| 733 |
+ |
|
| 734 |
+ |
nComponents = 1; |
| 735 |
+ |
component{ |
| 736 |
+ |
type = "Ar"; |
| 737 |
+ |
nMol = 108; |
| 738 |
+ |
} |
| 739 |
+ |
|
| 740 |
+ |
initialConfig = "./argon.init"; |
| 741 |
+ |
|
| 742 |
+ |
forceField = "LJ"; |
| 743 |
+ |
ensemble = "NVE"; // specify the simulation enesemble |
| 744 |
+ |
dt = 1.0; // the time step for integration |
| 745 |
+ |
runTime = 1e3; // the total simulation run time |
| 746 |
+ |
sampleTime = 100; // trajectory file frequency |
| 747 |
+ |
statusTime = 50; // statistics file frequency |
| 748 |
+ |
|
| 749 |
+ |
\end{lstlisting} |
| 750 |
+ |
|
| 751 |
|
Within the \texttt{.bass} file it is necessary to provide a complete |
| 752 |
|
description of the molecule before it is actually placed in the |
| 753 |
< |
simulation. The {\sc bass} syntax was originally developed with this goal in |
| 754 |
< |
mind, and allows for the specification of all the atoms in a molecular |
| 755 |
< |
prototype, as well as any bonds, bends, or torsions. These |
| 753 |
> |
simulation. The {\sc bass} syntax was originally developed with this |
| 754 |
> |
goal in mind, and allows for the specification of all the atoms in a |
| 755 |
> |
molecular prototype, as well as any bonds, bends, or torsions. These |
| 756 |
|
descriptions can become lengthy for complex molecules, and it would be |
| 757 |
< |
inconvenient to duplicate the simulation at the beginning of each {\sc bass} |
| 758 |
< |
script. Addressing this issue {\sc bass} allows for the inclusion of model |
| 759 |
< |
files at the top of a \texttt{.bass} file. These model files, denoted |
| 760 |
< |
with the \texttt{.mdl} extension, allow the user to describe a |
| 761 |
< |
molecular prototype once, then simply include it into each simulation |
| 762 |
< |
containing that molecule. |
| 757 |
> |
inconvenient to duplicate the simulation at the beginning of each {\sc |
| 758 |
> |
bass} script. Addressing this issue {\sc bass} allows for the |
| 759 |
> |
inclusion of model files at the top of a \texttt{.bass} file. These |
| 760 |
> |
model files, denoted with the \texttt{.mdl} extension, allow the user |
| 761 |
> |
to describe a molecular prototype once, then simply include it into |
| 762 |
> |
each simulation containing that molecule. Returning to the example in |
| 763 |
> |
Scheme~\ref{sch:bassExample}, the \texttt{.mdl} file's contents would |
| 764 |
> |
be Scheme~\ref{sch:mdlExample}, and the new \texttt{.bass} file would |
| 765 |
> |
become Scheme~\ref{sch:bassExPrime}. |
| 766 |
> |
|
| 767 |
> |
\begin{lstlisting}[float,caption={An example \texttt{.mdl} file.},label={sch:mdlExample}] |
| 768 |
> |
|
| 769 |
> |
molecule{ |
| 770 |
> |
name = "Ar"; |
| 771 |
> |
nAtoms = 1; |
| 772 |
> |
atom[0]{ |
| 773 |
> |
type="Ar"; |
| 774 |
> |
position( 0.0, 0.0, 0.0 ); |
| 775 |
> |
} |
| 776 |
> |
} |
| 777 |
> |
|
| 778 |
> |
\end{lstlisting} |
| 779 |
> |
|
| 780 |
> |
\begin{lstlisting}[float,caption={Revised {\sc bass} example.},label={sch:bassExPrime}] |
| 781 |
> |
|
| 782 |
> |
#include "argon.mdl" |
| 783 |
> |
|
| 784 |
> |
molecule{ |
| 785 |
> |
name = "Ar"; |
| 786 |
> |
nAtoms = 1; |
| 787 |
> |
atom[0]{ |
| 788 |
> |
type="Ar"; |
| 789 |
> |
position( 0.0, 0.0, 0.0 ); |
| 790 |
> |
} |
| 791 |
> |
} |
| 792 |
> |
|
| 793 |
> |
nComponents = 1; |
| 794 |
> |
component{ |
| 795 |
> |
type = "Ar"; |
| 796 |
> |
nMol = 108; |
| 797 |
> |
} |
| 798 |
> |
|
| 799 |
> |
initialConfig = "./argon.init"; |
| 800 |
> |
|
| 801 |
> |
forceField = "LJ"; |
| 802 |
> |
ensemble = "NVE"; |
| 803 |
> |
dt = 1.0; |
| 804 |
> |
runTime = 1e3; |
| 805 |
> |
sampleTime = 100; |
| 806 |
> |
statusTime = 50; |
| 807 |
|
|
| 808 |
+ |
\end{lstlisting} |
| 809 |
+ |
|
| 810 |
|
\subsection{\label{oopseSec:coordFiles}Coordinate Files} |
| 811 |
|
|
| 812 |
|
The standard format for storage of a systems coordinates is a modified |
| 813 |
|
xyz-file syntax, the exact details of which can be seen in |
| 814 |
< |
App.~\ref{appCoordFormat}. As all bonding and molecular information is |
| 815 |
< |
stored in the \texttt{.bass} and \texttt{.mdl} files, the coordinate |
| 816 |
< |
files are simply the complete set of coordinates for each atom at a |
| 817 |
< |
given simulation time. |
| 814 |
> |
Scheme~\ref{sch:dumpFormat}. As all bonding and molecular information |
| 815 |
> |
is stored in the \texttt{.bass} and \texttt{.mdl} files, the |
| 816 |
> |
coordinate files are simply the complete set of coordinates for each |
| 817 |
> |
atom at a given simulation time. One important note, although the |
| 818 |
> |
simulation propagates the complete rotation matrix, directional |
| 819 |
> |
entities are written out using quanternions, to save space in the |
| 820 |
> |
output files. |
| 821 |
|
|
| 822 |
< |
There are three major files used by {\sc oopse} written in the coordinate |
| 823 |
< |
format, they are as follows: the initialization file, the simulation |
| 824 |
< |
trajectory file, and the final coordinates of the simulation. The |
| 825 |
< |
initialization file is necessary for {\sc oopse} to start the simulation |
| 826 |
< |
with the proper coordinates. It is typically denoted with the |
| 827 |
< |
extension \texttt{.init}. The trajectory file is created at the |
| 828 |
< |
beginning of the simulation, and is used to store snapshots of the |
| 829 |
< |
simulation at regular intervals. The first frame is a duplication of |
| 830 |
< |
the \texttt{.init} file, and each subsequent frame is appended to the |
| 831 |
< |
file at an interval specified in the \texttt{.bass} file. The |
| 832 |
< |
trajectory file is given the extension \texttt{.dump}. The final |
| 833 |
< |
coordinate file is the end of run or \texttt{.eor} file. The |
| 822 |
> |
\begin{lstlisting}[float,caption={[The format of the coordinate files]Shows the format of the coordinate files. The fist line is the number of atoms. The second line begins with the time stamp followed by the three $\mathbf{H}$ column vectors. The next lines are the atomic coordinates for all atoms in the system. First is the name followed by position, velocity, quanternions, and lastly angular momentum.},label=sch:dumpFormat] |
| 823 |
> |
|
| 824 |
> |
nAtoms |
| 825 |
> |
time; Hxx Hyx Hzx; Hxy Hyy Hzy; Hxz Hyz Hzz; |
| 826 |
> |
Name1 x y z vx vy vz q0 q1 q2 q3 jx jy jz |
| 827 |
> |
Name2 x y z vx vy vz q0 q1 q2 q3 jx jy jz |
| 828 |
> |
etc... |
| 829 |
> |
|
| 830 |
> |
\end{lstlisting} |
| 831 |
> |
|
| 832 |
> |
|
| 833 |
> |
There are three major files used by {\sc oopse} written in the |
| 834 |
> |
coordinate format, they are as follows: the initialization file |
| 835 |
> |
(\texttt{.init}), the simulation trajectory file (\texttt{.dump}), and |
| 836 |
> |
the final coordinates of the simulation. The initialization file is |
| 837 |
> |
necessary for {\sc oopse} to start the simulation with the proper |
| 838 |
> |
coordinates, and is generated before the simulation run. The |
| 839 |
> |
trajectory file is created at the beginning of the simulation, and is |
| 840 |
> |
used to store snapshots of the simulation at regular intervals. The |
| 841 |
> |
first frame is a duplication of the |
| 842 |
> |
\texttt{.init} file, and each subsequent frame is appended to the file |
| 843 |
> |
at an interval specified in the \texttt{.bass} file with the |
| 844 |
> |
\texttt{sampleTime} flag. The final coordinate file is the end of run file. The |
| 845 |
|
\texttt{.eor} file stores the final configuration of the system for a |
| 846 |
|
given simulation. The file is updated at the same time as the |
| 847 |
< |
\texttt{.dump} file. However, it only contains the most recent |
| 847 |
> |
\texttt{.dump} file, however, it only contains the most recent |
| 848 |
|
frame. In this way, an \texttt{.eor} file may be used as the |
| 849 |
< |
initialization file to a second simulation in order to continue or |
| 850 |
< |
recover the previous simulation. |
| 849 |
> |
initialization file to a second simulation in order to continue a |
| 850 |
> |
simulation or recover one from a processor that has crashed during the |
| 851 |
> |
course of the run. |
| 852 |
|
|
| 853 |
|
\subsection{\label{oopseSec:initCoords}Generation of Initial Coordinates} |
| 854 |
|
|
| 855 |
< |
As was stated in Sec.~\ref{oopseSec:coordFiles}, an initialization file |
| 856 |
< |
is needed to provide the starting coordinates for a simulation. The |
| 857 |
< |
{\sc oopse} package provides a program called \texttt{sysBuilder} to aid in |
| 858 |
< |
the creation of the \texttt{.init} file. \texttt{sysBuilder} is {\sc bass} |
| 859 |
< |
aware, and will recognize arguments and parameters in the |
| 860 |
< |
\texttt{.bass} file that would otherwise be ignored by the |
| 861 |
< |
simulation. The program itself is under continual development, and is |
| 782 |
< |
offered here as a helper tool only. |
| 855 |
> |
As was stated in Sec.~\ref{oopseSec:coordFiles}, an initialization |
| 856 |
> |
file is needed to provide the starting coordinates for a |
| 857 |
> |
simulation. The {\sc oopse} package provides a program called |
| 858 |
> |
\texttt{sysBuilder} to aid in the creation of the \texttt{.init} |
| 859 |
> |
file. \texttt{sysBuilder} uses {\sc bass}, and will recognize |
| 860 |
> |
arguments and parameters in the \texttt{.bass} file that would |
| 861 |
> |
otherwise be ignored by the simulation. |
| 862 |
|
|
| 863 |
|
\subsection{The Statistics File} |
| 864 |
|
|
| 865 |
< |
The last output file generated by {\sc oopse} is the statistics file. This |
| 866 |
< |
file records such statistical quantities as the instantaneous |
| 867 |
< |
temperature, volume, pressure, etc. It is written out with the |
| 868 |
< |
frequency specified in the \texttt{.bass} file. The file allows the |
| 869 |
< |
user to observe the system variables as a function of simulation time |
| 870 |
< |
while the simulation is in progress. One useful function the |
| 871 |
< |
statistics file serves is to monitor the conserved quantity of a given |
| 872 |
< |
simulation ensemble, this allows the user to observe the stability of |
| 873 |
< |
the integrator. The statistics file is denoted with the \texttt{.stat} |
| 874 |
< |
file extension. |
| 865 |
> |
The last output file generated by {\sc oopse} is the statistics |
| 866 |
> |
file. This file records such statistical quantities as the |
| 867 |
> |
instantaneous temperature, volume, pressure, etc. It is written out |
| 868 |
> |
with the frequency specified in the \texttt{.bass} file with the |
| 869 |
> |
\texttt{statusTime} keyword. The file allows the user to observe the |
| 870 |
> |
system variables as a function of simulation time while the simulation |
| 871 |
> |
is in progress. One useful function the statistics file serves is to |
| 872 |
> |
monitor the conserved quantity of a given simulation ensemble, this |
| 873 |
> |
allows the user to observe the stability of the integrator. The |
| 874 |
> |
statistics file is denoted with the \texttt{.stat} file extension. |
| 875 |
|
|
| 876 |
|
\section{\label{oopseSec:mechanics}Mechanics} |
| 877 |
|
|
| 879 |
|
|
| 880 |
|
Integration of the equations of motion was carried out using the |
| 881 |
|
symplectic splitting method proposed by Dullweber \emph{et |
| 882 |
< |
al.}.\cite{Dullweber1997} The reason for this integrator selection |
| 883 |
< |
deals with poor energy conservation of rigid body systems using |
| 884 |
< |
quaternions. While quaternions work well for orientational motion in |
| 885 |
< |
alternate ensembles, the microcanonical ensemble has a constant energy |
| 886 |
< |
requirement that is quite sensitive to errors in the equations of |
| 887 |
< |
motion. The original implementation of this code utilized quaternions |
| 888 |
< |
for rotational motion propagation; however, a detailed investigation |
| 889 |
< |
showed that they resulted in a steady drift in the total energy, |
| 890 |
< |
something that has been observed by others.\cite{Laird97} |
| 882 |
> |
al.}.\cite{Dullweber1997} The reason for the selection of this |
| 883 |
> |
integrator, is the poor energy conservation of rigid body systems |
| 884 |
> |
using quaternion dynamics. While quaternions work well for |
| 885 |
> |
orientational motion in alternate ensembles, the microcanonical |
| 886 |
> |
ensemble has a constant energy requirement that is quite sensitive to |
| 887 |
> |
errors in the equations of motion. The original implementation of {\sc |
| 888 |
> |
oopse} utilized quaternions for rotational motion propagation; |
| 889 |
> |
however, a detailed investigation showed that they resulted in a |
| 890 |
> |
steady drift in the total energy, something that has been observed by |
| 891 |
> |
others.\cite{Laird97} |
| 892 |
|
|
| 893 |
|
The key difference in the integration method proposed by Dullweber |
| 894 |
< |
\emph{et al.} is that the entire rotation matrix is propagated from |
| 895 |
< |
one time step to the next. In the past, this would not have been as |
| 896 |
< |
feasible a option, being that the rotation matrix for a single body is |
| 897 |
< |
nine elements long as opposed to 3 or 4 elements for Euler angles and |
| 898 |
< |
quaternions respectively. System memory has become much less of an |
| 899 |
< |
issue in recent times, and this has resulted in substantial benefits |
| 900 |
< |
in energy conservation. There is still the issue of 5 or 6 additional |
| 821 |
< |
elements for describing the orientation of each particle, which will |
| 822 |
< |
increase dump files substantially. Simply translating the rotation |
| 823 |
< |
matrix into its component Euler angles or quaternions for storage |
| 824 |
< |
purposes relieves this burden. |
| 894 |
> |
\emph{et al}.~({\sc dlm}) is that the entire rotation matrix is propagated from |
| 895 |
> |
one time step to the next. In the past, this would not have been a |
| 896 |
> |
feasible option, since the rotation matrix for a single body is nine |
| 897 |
> |
elements long as opposed to three or four elements for Euler angles |
| 898 |
> |
and quaternions respectively. System memory has become much less of an |
| 899 |
> |
issue in recent times, and the {\sc dlm} method has used memory in |
| 900 |
> |
exchange for substantial benefits in energy conservation. |
| 901 |
|
|
| 902 |
< |
The symplectic splitting method allows for Verlet style integration of |
| 903 |
< |
both linear and angular motion of rigid bodies. In the integration |
| 904 |
< |
method, the orientational propagation involves a sequence of matrix |
| 902 |
> |
The {\sc dlm} method allows for Verlet style integration of both |
| 903 |
> |
linear and angular motion of rigid bodies. In the integration method, |
| 904 |
> |
the orientational propagation involves a sequence of matrix |
| 905 |
|
evaluations to update the rotation matrix.\cite{Dullweber1997} These |
| 906 |
< |
matrix rotations end up being more costly computationally than the |
| 907 |
< |
simpler arithmetic quaternion propagation. With the same time step, a |
| 908 |
< |
1000 SSD particle simulation shows an average 7\% increase in |
| 909 |
< |
computation time using the symplectic step method in place of |
| 910 |
< |
quaternions. This cost is more than justified when comparing the |
| 911 |
< |
energy conservation of the two methods as illustrated in figure |
| 836 |
< |
\ref{timestep}. |
| 906 |
> |
matrix rotations are more costly computationally than the simpler |
| 907 |
> |
arithmetic quaternion propagation. With the same time step, a 1000 SSD |
| 908 |
> |
particle simulation shows an average 7\% increase in computation time |
| 909 |
> |
using the {\sc dlm} method in place of quaternions. This cost is more |
| 910 |
> |
than justified when comparing the energy conservation of the two |
| 911 |
> |
methods as illustrated in Fig.~\ref{timestep}. |
| 912 |
|
|
| 913 |
|
\begin{figure} |
| 914 |
|
\centering |
| 915 |
|
\includegraphics[width=\linewidth]{timeStep.eps} |
| 916 |
< |
\caption{Energy conservation using quaternion based integration versus |
| 917 |
< |
the symplectic step method proposed by Dullweber \emph{et al.} with |
| 916 |
> |
\caption[Energy conservation for quaternion versus {\sc dlm} dynamics]{Energy conservation using quaternion based integration versus |
| 917 |
> |
the {\sc dlm} method with |
| 918 |
|
increasing time step. For each time step, the dotted line is total |
| 919 |
< |
energy using the symplectic step integrator, and the solid line comes |
| 919 |
> |
energy using the {\sc dlm} integrator, and the solid line comes |
| 920 |
|
from the quaternion integrator. The larger time step plots are shifted |
| 921 |
|
up from the true energy baseline for clarity.} |
| 922 |
|
\label{timestep} |
| 923 |
|
\end{figure} |
| 924 |
|
|
| 925 |
< |
In figure \ref{timestep}, the resulting energy drift at various time |
| 926 |
< |
steps for both the symplectic step and quaternion integration schemes |
| 925 |
> |
In Fig.~\ref{timestep}, the resulting energy drift at various time |
| 926 |
> |
steps for both the {\sc dlm} and quaternion integration schemes |
| 927 |
|
is compared. All of the 1000 SSD particle simulations started with the |
| 928 |
|
same configuration, and the only difference was the method for |
| 929 |
|
handling rotational motion. At time steps of 0.1 and 0.5 fs, both |
| 930 |
|
methods for propagating particle rotation conserve energy fairly well, |
| 931 |
|
with the quaternion method showing a slight energy drift over time in |
| 932 |
|
the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the |
| 933 |
< |
energy conservation benefits of the symplectic step method are clearly |
| 933 |
> |
energy conservation benefits of the {\sc dlm} method are clearly |
| 934 |
|
demonstrated. Thus, while maintaining the same degree of energy |
| 935 |
|
conservation, one can take considerably longer time steps, leading to |
| 936 |
|
an overall reduction in computation time. |
| 939 |
|
time steps up to three femtoseconds. A slight energy drift on the |
| 940 |
|
order of 0.012 kcal/mol per nanosecond was observed at a time step of |
| 941 |
|
four femtoseconds, and as expected, this drift increases dramatically |
| 942 |
< |
with increasing time step. To insure accuracy in the constant energy |
| 868 |
< |
simulations, time steps were set at 2 fs and kept at this value for |
| 869 |
< |
constant pressure simulations as well. |
| 942 |
> |
with increasing time step. |
| 943 |
|
|
| 944 |
|
|
| 945 |
|
\subsection{\label{sec:extended}Extended Systems for other Ensembles} |
| 948 |
|
{\sc oopse} implements a |
| 949 |
|
|
| 950 |
|
|
| 951 |
< |
\subsubsection{\label{oopseSec:noseHooverThermo}Nose-Hoover Thermostatting} |
| 951 |
> |
\subsection{\label{oopseSec:noseHooverThermo}Nose-Hoover Thermostatting} |
| 952 |
|
|
| 953 |
|
To mimic the effects of being in a constant temperature ({\sc nvt}) |
| 954 |
|
ensemble, {\sc oopse} uses the Nose-Hoover extended system |
| 974 |
|
some subtlety in choosing values for $\tau_{T}$, and it is usually set |
| 975 |
|
to values of a few ps. Within a {\sc bass} file, $\tau_{T}$ could be |
| 976 |
|
set to 1 ps using the {\tt tauThermostat = 1000; } command. |
| 977 |
+ |
|
| 978 |
+ |
\subsection{\label{oopseSec:rattle}The {\sc rattle} Method for Bond |
| 979 |
+ |
Constraints} |
| 980 |
+ |
|
| 981 |
+ |
In order to satisfy the constraints of fixed bond lengths within {\sc |
| 982 |
+ |
oopse}, we have implemented the {\sc rattle} algorithm of |
| 983 |
+ |
Andersen.\cite{andersen83} The algorithm is a velocity verlet |
| 984 |
+ |
formulation of the {\sc shake} method\cite{ryckaert77} of iteratively |
| 985 |
+ |
solving the Lagrange multipliers of constraint. The system of lagrange |
| 986 |
+ |
multipliers allows one to reformulate the equations of motion with |
| 987 |
+ |
explicit constraint forces on the equations of |
| 988 |
+ |
motion.\cite{fowles99:lagrange} |
| 989 |
+ |
|
| 990 |
+ |
Consider a system described by qoordinates $q_1$ and $q_2$ subject to an |
| 991 |
+ |
equation of constraint: |
| 992 |
+ |
\begin{equation} |
| 993 |
+ |
\sigma(q_1, q_2,t) = 0 |
| 994 |
+ |
\label{oopseEq:lm1} |
| 995 |
+ |
\end{equation} |
| 996 |
+ |
The Lagrange formulation of the equations of motion can be written: |
| 997 |
+ |
\begin{equation} |
| 998 |
+ |
\delta\int_{t_1}^{t_2}L\, dt = |
| 999 |
+ |
\int_{t_1}^{t_2} \sum_i \biggl [ \frac{\partial L}{\partial q_i} |
| 1000 |
+ |
- \frac{d}{dt}\biggl(\frac{\partial L}{\partial \dot{q}_i} |
| 1001 |
+ |
\biggr ) \biggr] \delta q_i \, dt = 0 |
| 1002 |
+ |
\label{oopseEq:lm2} |
| 1003 |
+ |
\end{equation} |
| 1004 |
+ |
Here, $\delta q_i$ is not independent for each $q$, as $q_1$ and $q_2$ |
| 1005 |
+ |
are linked by $\sigma$. However, $\sigma$ is fixed at any given |
| 1006 |
+ |
instant of time, giving: |
| 1007 |
+ |
\begin{align} |
| 1008 |
+ |
\delta\sigma &= \biggl( \frac{\partial\sigma}{\partial q_1} \delta q_1 % |
| 1009 |
+ |
+ \frac{\partial\sigma}{\partial q_2} \delta q_2 \biggr) = 0 \\ |
| 1010 |
+ |
% |
| 1011 |
+ |
\frac{\partial\sigma}{\partial q_1} \delta q_1 &= % |
| 1012 |
+ |
- \frac{\partial\sigma}{\partial q_2} \delta q_2 \\ |
| 1013 |
+ |
% |
| 1014 |
+ |
\delta q_2 &= - \biggl(\frac{\partial\sigma}{\partial q_1} \bigg / % |
| 1015 |
+ |
\frac{\partial\sigma}{\partial q_2} \biggr) \delta q_1 |
| 1016 |
+ |
\end{align} |
| 1017 |
+ |
Substituted back into Eq.~\ref{oopseEq:lm2}, |
| 1018 |
+ |
\begin{equation} |
| 1019 |
+ |
\int_{t_1}^{t_2}\biggl [ \biggl(\frac{\partial L}{\partial q_1} |
| 1020 |
+ |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1021 |
+ |
\biggr) |
| 1022 |
+ |
- \biggl( \frac{\partial L}{\partial q_1} |
| 1023 |
+ |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1024 |
+ |
\biggr) \biggl(\frac{\partial\sigma}{\partial q_1} \bigg / % |
| 1025 |
+ |
\frac{\partial\sigma}{\partial q_2} \biggr)\biggr] \delta q_1 \, dt = 0 |
| 1026 |
+ |
\label{oopseEq:lm3} |
| 1027 |
+ |
\end{equation} |
| 1028 |
+ |
Leading to, |
| 1029 |
+ |
\begin{equation} |
| 1030 |
+ |
\frac{\biggl(\frac{\partial L}{\partial q_1} |
| 1031 |
+ |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1032 |
+ |
\biggr)}{\frac{\partial\sigma}{\partial q_1}} = |
| 1033 |
+ |
\frac{\biggl(\frac{\partial L}{\partial q_2} |
| 1034 |
+ |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_2} |
| 1035 |
+ |
\biggr)}{\frac{\partial\sigma}{\partial q_2}} |
| 1036 |
+ |
\label{oopseEq:lm4} |
| 1037 |
+ |
\end{equation} |
| 1038 |
+ |
This relation can only be statisfied, if both are equal to a single |
| 1039 |
+ |
function $-\lambda(t)$, |
| 1040 |
+ |
\begin{align} |
| 1041 |
+ |
\frac{\biggl(\frac{\partial L}{\partial q_1} |
| 1042 |
+ |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1043 |
+ |
\biggr)}{\frac{\partial\sigma}{\partial q_1}} &= -\lambda(t) \\ |
| 1044 |
+ |
% |
| 1045 |
+ |
\frac{\partial L}{\partial q_1} |
| 1046 |
+ |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} &= |
| 1047 |
+ |
-\lambda(t)\,\frac{\partial\sigma}{\partial q_1} \\ |
| 1048 |
+ |
% |
| 1049 |
+ |
\frac{\partial L}{\partial q_1} |
| 1050 |
+ |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1051 |
+ |
+ \mathcal{G}_i &= 0 |
| 1052 |
+ |
\end{align} |
| 1053 |
+ |
Where $\mathcal{G}_i$, the force of constraint on $i$, is: |
| 1054 |
+ |
\begin{equation} |
| 1055 |
+ |
\mathcal{G}_i = \lambda(t)\,\frac{\partial\sigma}{\partial q_1} |
| 1056 |
+ |
\label{oopseEq:lm5} |
| 1057 |
+ |
\end{equation} |
| 1058 |
|
|
| 1059 |
+ |
In a simulation, this would involve the solution of a set of $(m + n)$ |
| 1060 |
+ |
number of equations. Where $m$ is the number of constraints, and $n$ |
| 1061 |
+ |
is the number of constrained coordinates. In practice, this is not |
| 1062 |
+ |
done, as the matrix inversion neccassary to solve the system of |
| 1063 |
+ |
equations would be very time consuming to solve. Additionally, the |
| 1064 |
+ |
numerical error in the solution of the set of $\lambda$'s would be |
| 1065 |
+ |
compounded by the error inherent in propagating by the Velocity Verlet |
| 1066 |
+ |
algorithm ($\Delta t^4$). The verlet propagation error is negligible |
| 1067 |
+ |
in an unconstrained system, as one is interested in the statisitics of |
| 1068 |
+ |
the run, and not that the run be numerically exact to the ``true'' |
| 1069 |
+ |
integration. This relates back to the ergodic hypothesis that a time |
| 1070 |
+ |
integral of a valid trajectory will still give the correct enesemble |
| 1071 |
+ |
average. However, in the case of constraints, if the equations of |
| 1072 |
+ |
motion leave the ``true'' trajectory, they are departing from the |
| 1073 |
+ |
constrained surface. The method that is used, is to iteratively solve |
| 1074 |
+ |
for $\lambda(t)$ at each time step. |
| 1075 |
+ |
|
| 1076 |
+ |
In {\sc rattle} the equations of motion are modified subject to the |
| 1077 |
+ |
following two constraints: |
| 1078 |
+ |
\begin{align} |
| 1079 |
+ |
\sigma_{ij}[\mathbf{r}(t)] \equiv |
| 1080 |
+ |
[ \mathbf{r}_i(t) - \mathbf{r}_j(t)]^2 - d_{ij}^2 &= 0 % |
| 1081 |
+ |
\label{oopseEq:c1} \\ |
| 1082 |
+ |
% |
| 1083 |
+ |
[\mathbf{\dot{r}}_i(t) - \mathbf{\dot{r}}_j(t)] \cdot |
| 1084 |
+ |
[\mathbf{r}_i(t) - \mathbf{r}_j(t)] &= 0 \label{oopseEq:c2} |
| 1085 |
+ |
\end{align} |
| 1086 |
+ |
Eq.~\ref{oopseEq:c1} is the set of bond constraints, where $d_{ij}$ is |
| 1087 |
+ |
the constrained distance between atom $i$ and |
| 1088 |
+ |
$j$. Eq.~\ref{oopseEq:c2} constrains the velocities of $i$ and $j$ to |
| 1089 |
+ |
be perpindicular to the bond vector, so that the bond can neither grow |
| 1090 |
+ |
nor shrink. The constrained dynamics equations become: |
| 1091 |
+ |
\begin{equation} |
| 1092 |
+ |
m_i \mathbf{\ddot{r}}_i = \mathbf{F}_i + \mathbf{\mathcal{G}}_i |
| 1093 |
+ |
\label{oopseEq:r1} |
| 1094 |
+ |
\end{equation} |
| 1095 |
+ |
Where, |
| 1096 |
+ |
\begin{equation} |
| 1097 |
+ |
\mathbf{\mathcal{G}}_i = - \sum_j \lambda_{ij}(t)\,\nabla \sigma_{ij} |
| 1098 |
+ |
\label{oopseEq:r2} |
| 1099 |
+ |
\end{equation} |
| 1100 |
+ |
|
| 1101 |
+ |
In Velocity Verlet, if $\Delta t = h$, the propagation can be written: |
| 1102 |
+ |
\begin{align} |
| 1103 |
+ |
\mathbf{r}_i(t+h) &= |
| 1104 |
+ |
\mathbf{r}_i(t) + h\mathbf{\dot{r}}(t) + |
| 1105 |
+ |
\frac{h^2}{2m_i}\,\Bigl[ \mathbf{F}_i(t) + |
| 1106 |
+ |
\mathbf{\mathcal{G}}_{Ri}(t) \Bigr] \label{oopseEq:vv1} \\ |
| 1107 |
+ |
% |
| 1108 |
+ |
\mathbf{\dot{r}}_i(t+h) &= |
| 1109 |
+ |
\mathbf{\dot{r}}_i(t) + \frac{h}{2m_i} |
| 1110 |
+ |
\Bigl[ \mathbf{F}_i(t) + \mathbf{\mathcal{G}}_{Ri}(t) + |
| 1111 |
+ |
\mathbf{F}_i(t+h) + \mathbf{\mathcal{G}}_{Vi}(t+h) \Bigr] % |
| 1112 |
+ |
\label{oopseEq:vv2} |
| 1113 |
+ |
\end{align} |
| 1114 |
|
|
| 1115 |
+ |
|
| 1116 |
+ |
|
| 1117 |
|
\subsection{\label{oopseSec:zcons}Z-Constraint Method} |
| 1118 |
|
|
| 1119 |
< |
Based on fluctuation-dissipation theorem,\bigskip\ force auto-correlation |
| 1119 |
> |
Based on fluctuation-dissipation theorem, a force auto-correlation |
| 1120 |
|
method was developed to investigate the dynamics of ions inside the ion |
| 1121 |
|
channels.\cite{Roux91} Time-dependent friction coefficient can be calculated |
| 1122 |
|
from the deviation of the instantaneous force from its mean force. |
