ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mattDisertation/Introduction.tex
(Generate patch)

Comparing trunk/mattDisertation/Introduction.tex (file contents):
Revision 1009 by mmeineke, Tue Feb 3 17:54:39 2004 UTC vs.
Revision 1012 by mmeineke, Tue Feb 3 21:14:39 2004 UTC

# Line 568 | Line 568 | q(t+\Delta t) =
568   \end{equation}
569   Or equivalently,
570   \begin{equation}
571 < q(t+\Delta t) =
572 <        2q(t) - q(t-\Delta t) + \frac{F(t)}{m}\Delta t^2 +
573 <        \mathcal{O}(\Delta t^4)
571 > q(t+\Delta t) \approx
572 >        2q(t) - q(t-\Delta t) + \frac{F(t)}{m}\Delta t^2
573   \label{introEq:verletFinal}
574   \end{equation}
575   Which contains an error in the estimate of the new positions on the
# Line 578 | Line 577 | with a velocity reformulation of the Verlet method.\ci
577  
578   In practice, however, the simulations in this research were integrated
579   with a velocity reformulation of the Verlet method.\cite{allen87:csl}
580 < \begin{equation}
581 < q(t+\Delta t)= q(t) + v(t)\Delta t + \frac{F(t)}{2m}\Delta t^2
582 < \label{introEq:MDvelVerletPos}
583 < \end{equation}
584 < \begin{equation}
586 < v(t+\Delta t) = v(t) + \frac{\Delta t}{2m}[F(t) + F(t+\Delta t)]
580 > \begin{align}
581 > q(t+\Delta t) &= q(t) + v(t)\Delta t + \frac{F(t)}{2m}\Delta t^2 %
582 > \label{introEq:MDvelVerletPos} \\%
583 > %
584 > v(t+\Delta t) &= v(t) + \frac{\Delta t}{2m}[F(t) + F(t+\Delta t)] %
585   \label{introEq:MDvelVerletVel}
586 < \end{equation}
586 > \end{align}
587   The original Verlet algorithm can be regained by substituting the
588   velocity back into Eq.~\ref{introEq:MDvelVerletPos}.  The Verlet
589   formulations are chosen in this research because the algorithms have
# Line 607 | Line 605 | hypothesis (Sec.~\ref{introSec:StatThermo}), it is kno
605   reversible.  The fact that it shadows the true Hamiltonian in phase
606   space is acceptable in actual simulations as one is interested in the
607   ensemble average of the observable being measured.  From the ergodic
608 < hypothesis (Sec.~\ref{introSec:StatThermo}), it is known that the time
608 > hypothesis (Sec.~\ref{introSec:statThermo}), it is known that the time
609   average will match the ensemble average, therefore two similar
610   trajectories in phase space should give matching statistical averages.
611  
612   \subsection{\label{introSec:MDfurther}Further Considerations}
613 +
614   In the simulations presented in this research, a few additional
615   parameters are needed to describe the motions.  The simulations
616 < involving water and phospholipids in Ch.~\ref{chaptLipids} are
616 > involving water and phospholipids in Ch.~\ref{chapt:lipid} are
617   required to integrate the equations of motions for dipoles on atoms.
618   This involves an additional three parameters be specified for each
619   dipole atom: $\phi$, $\theta$, and $\psi$.  These three angles are
620   taken to be the Euler angles, where $\phi$ is a rotation about the
621   $z$-axis, and $\theta$ is a rotation about the new $x$-axis, and
622   $\psi$ is a final rotation about the new $z$-axis (see
623 < Fig.~\ref{introFig:euleerAngles}).  This sequence of rotations can be
624 < accumulated into a single $3 \times 3$ matrix $\mathbf{A}$
623 > Fig.~\ref{introFig:eulerAngles}).  This sequence of rotations can be
624 > accumulated into a single $3 \times 3$ matrix, $\mathbf{A}$,
625   defined as follows:
626   \begin{equation}
627 < eq here
627 > \mathbf{A} =
628 > \begin{bmatrix}
629 >        \cos\phi\cos\psi-\sin\phi\cos\theta\sin\psi &%
630 >        \sin\phi\cos\psi+\cos\phi\cos\theta\sin\psi &%
631 >        \sin\theta\sin\psi \\%
632 >        %
633 >        -\cos\phi\sin\psi-\sin\phi\cos\theta\cos\psi &%
634 >        -\sin\phi\sin\psi+\cos\phi\cos\theta\cos\psi &%
635 >        \sin\theta\cos\psi \\%
636 >        %
637 >        \sin\phi\sin\theta &%
638 >        -\cos\phi\sin\theta &%
639 >        \cos\theta
640 > \end{bmatrix}
641   \label{introEq:EulerRotMat}
642   \end{equation}
643  
644 < The equations of motion for Euler angles can be written down as
645 < \cite{allen87:csl}
646 < \begin{equation}
647 < eq here
648 < \label{introEq:MDeuleeerPsi}
649 < \end{equation}
644 > The equations of motion for Euler angles can be written down
645 > as\cite{allen87:csl}
646 > \begin{align}
647 > \dot{\phi} &= -\omega^s_x \frac{\sin\phi\cos\theta}{\sin\theta} +
648 >        \omega^s_y \frac{\cos\phi\cos\theta}{\sin\theta} +
649 >        \omega^s_z
650 > \label{introEq:MDeulerPhi} \\%
651 > %
652 > \dot{\theta} &= \omega^s_x \cos\phi + \omega^s_y \sin\phi
653 > \label{introEq:MDeulerTheta} \\%
654 > %
655 > \dot{\psi} &= \omega^s_x \frac{\sin\phi}{\sin\theta} -
656 >        \omega^s_y \frac{\cos\phi}{\sin\theta}
657 > \label{introEq:MDeulerPsi}
658 > \end{align}
659   Where $\omega^s_i$ is the angular velocity in the lab space frame
660   along Cartesian coordinate $i$.  However, a difficulty arises when
661   attempting to integrate Eq.~\ref{introEq:MDeulerPhi} and
662   Eq.~\ref{introEq:MDeulerPsi}. The $\frac{1}{\sin \theta}$ present in
663   both equations means there is a non-physical instability present when
664 < $\theta$ is 0 or $\pi$.
665 <
666 < To correct for this, the simulations integrate the rotation matrix,
667 < $\mathbf{A}$, directly, thus avoiding the instability.
647 < This method was proposed by Dullwebber
648 < \emph{et. al.}\cite{Dullwebber:1997}, and is presented in
664 > $\theta$ is 0 or $\pi$. To correct for this, the simulations integrate
665 > the rotation matrix, $\mathbf{A}$, directly, thus avoiding the
666 > instability.  This method was proposed by Dullweber
667 > \emph{et. al.}\cite{Dullweber1997}, and is presented in
668   Sec.~\ref{introSec:MDsymplecticRot}.
669  
670 < \subsubsection{\label{introSec:MDliouville}Liouville Propagator}
670 > \subsection{\label{introSec:MDliouville}Liouville Propagator}
671  
672   Before discussing the integration of the rotation matrix, it is
673   necessary to understand the construction of a ``good'' integration
674   scheme.  It has been previously
675 < discussed(Sec.~\ref{introSec:MDintegrate}) how it is desirable for an
675 > discussed(Sec.~\ref{introSec:mdIntegrate}) how it is desirable for an
676   integrator to be symplectic, or time reversible.  The following is an
677   outline of the Trotter factorization of the Liouville Propagator as a
678 < scheme for generating symplectic integrators. \cite{Tuckerman:1992}
678 > scheme for generating symplectic integrators.\cite{Tuckerman92}
679  
680   For a system with $f$ degrees of freedom the Liouville operator can be
681   defined as,
682   \begin{equation}
683 < eq here
683 > iL=\sum^f_{j=1} \biggl [\dot{q}_j\frac{\partial}{\partial q_j} +
684 >        F_j\frac{\partial}{\partial p_j} \biggr ]
685   \label{introEq:LiouvilleOperator}
686   \end{equation}
687 < Here, $r_j$ and $p_j$ are the position and conjugate momenta of a
688 < degree of freedom, and $f_j$ is the force on that degree of freedom.
687 > Here, $q_j$ and $p_j$ are the position and conjugate momenta of a
688 > degree of freedom, and $F_j$ is the force on that degree of freedom.
689   $\Gamma$ is defined as the set of all positions and conjugate momenta,
690 < $\{r_j,p_j\}$, and the propagator, $U(t)$, is defined
690 > $\{q_j,p_j\}$, and the propagator, $U(t)$, is defined
691   \begin {equation}
692 < eq here
692 > U(t) = e^{iLt}
693   \label{introEq:Lpropagator}
694   \end{equation}
695   This allows the specification of $\Gamma$ at any time $t$ as
696   \begin{equation}
697 < eq here
697 > \Gamma(t) = U(t)\Gamma(0)
698   \label{introEq:Lp2}
699   \end{equation}
700   It is important to note, $U(t)$ is a unitary operator meaning
# Line 685 | Line 705 | Trotter theorem to yield
705  
706   Decomposing $L$ into two parts, $iL_1$ and $iL_2$, one can use the
707   Trotter theorem to yield
708 < \begin{equation}
709 < eq here
710 < \label{introEq:Lp4}
711 < \end{equation}
712 < Where $\Delta t = \frac{t}{P}$.
708 > \begin{align}
709 > e^{iLt} &= e^{i(L_1 + L_2)t} \notag \\%
710 > %
711 >        &= \biggl [ e^{i(L_1 +L_2)\frac{t}{P}} \biggr]^P \notag \\%
712 > %
713 >        &= \biggl [ e^{iL_1\frac{\Delta t}{2}}\, e^{iL_2\Delta t}\,
714 >        e^{iL_1\frac{\Delta t}{2}} \biggr ]^P +
715 >        \mathcal{O}\biggl (\frac{t^3}{P^2} \biggr ) \label{introEq:Lp4}
716 > \end{align}
717 > Where $\Delta t = t/P$.
718   With this, a discrete time operator $G(\Delta t)$ can be defined:
719 < \begin{equation}
720 < eq here
719 > \begin{align}
720 > G(\Delta t) &= e^{iL_1\frac{\Delta t}{2}}\, e^{iL_2\Delta t}\,
721 >        e^{iL_1\frac{\Delta t}{2}} \notag \\%
722 > %
723 >        &= U_1 \biggl ( \frac{\Delta t}{2} \biggr )\, U_2 ( \Delta t )\,
724 >        U_1 \biggl ( \frac{\Delta t}{2} \biggr )
725   \label{introEq:Lp5}
726 < \end{equation}
727 < Because $U_1(t)$ and $U_2(t)$ are unitary, $G|\Delta t)$ is also
726 > \end{align}
727 > Because $U_1(t)$ and $U_2(t)$ are unitary, $G(\Delta t)$ is also
728   unitary.  Meaning an integrator based on this factorization will be
729   reversible in time.
730  
731   As an example, consider the following decomposition of $L$:
732 + \begin{align}
733 + iL_1 &= \dot{q}\frac{\partial}{\partial q}%
734 + \label{introEq:Lp6a} \\%
735 + %
736 + iL_2 &= F(q)\frac{\partial}{\partial p}%
737 + \label{introEq:Lp6b}
738 + \end{align}
739 + This leads to propagator $G( \Delta t )$ as,
740   \begin{equation}
741 < eq here
742 < \label{introEq:Lp6}
741 > G(\Delta t) =  e^{\frac{\Delta t}{2} F(q)\frac{\partial}{\partial p}} \,
742 >        e^{\Delta t\,\dot{q}\frac{\partial}{\partial q}} \,
743 >        e^{\frac{\Delta t}{2} F(q)\frac{\partial}{\partial p}}
744 > \label{introEq:Lp7}
745   \end{equation}
746 < Operating $G(\Delta t)$ on $\Gamma)0)$, and utilizing the operator property
746 > Operating $G(\Delta t)$ on $\Gamma(0)$, and utilizing the operator property
747   \begin{equation}
748 < eq here
748 > e^{c\frac{\partial}{\partial x}}\, f(x) = f(x+c)
749   \label{introEq:Lp8}
750   \end{equation}
751 < Where $c$ is independent of $q$.  One obtains the following:  
752 < \begin{equation}
753 < eq here
754 < \label{introEq:Lp8}
755 < \end{equation}
751 > Where $c$ is independent of $x$.  One obtains the following:  
752 > \begin{align}
753 > \dot{q}\biggl (\frac{\Delta t}{2}\biggr ) &=
754 >        \dot{q}(0) + \frac{\Delta t}{2m}\, F[q(0)] \label{introEq:Lp9a}\\%
755 > %
756 > q(\Delta t) &= q(0) + \Delta t\, \dot{q}\biggl (\frac{\Delta t}{2}\biggr )%
757 >        \label{introEq:Lp9b}\\%
758 > %
759 > \dot{q}(\Delta t) &= \dot{q}\biggl (\frac{\Delta t}{2}\biggr ) +
760 >        \frac{\Delta t}{2m}\, F[q(0)] \label{introEq:Lp9c}
761 > \end{align}
762   Or written another way,
763 < \begin{equation}
764 < eq here
765 < \label{intorEq:Lp9}
766 < \end{equation}
763 > \begin{align}
764 > q(t+\Delta t) &= q(0) + \dot{q}(0)\Delta t +
765 >        \frac{F[q(0)]}{m}\frac{\Delta t^2}{2} %
766 > \label{introEq:Lp10a} \\%
767 > %
768 > \dot{q}(\Delta t) &= \dot{q}(0) + \frac{\Delta t}{2m}
769 >        \biggl [F[q(0)] + F[q(\Delta t)] \biggr] %
770 > \label{introEq:Lp10b}
771 > \end{align}
772   This is the velocity Verlet formulation presented in
773 < Sec.~\ref{introSec:MDintegrate}.  Because this integration scheme is
773 > Sec.~\ref{introSec:mdIntegrate}.  Because this integration scheme is
774   comprised of unitary propagators, it is symplectic, and therefore area
775   preserving in phase space.  From the preceding factorization, one can
776   see that the integration of the equations of motion would follow:
# Line 734 | Line 784 | see that the integration of the equations of motion wo
784   \item Repeat from step 1 with the new position, velocities, and forces assuming the roles of the initial values.
785   \end{enumerate}
786  
787 < \subsubsection{\label{introSec:MDsymplecticRot} Symplectic Propagation of the Rotation Matrix}
787 > \subsection{\label{introSec:MDsymplecticRot} Symplectic Propagation of the Rotation Matrix}
788  
789   Based on the factorization from the previous section,
790 < Dullweber\emph{et al.}\cite{Dullweber:1997}~ proposed a scheme for the
790 > Dullweber\emph{et al}.\cite{Dullweber1997}~ proposed a scheme for the
791   symplectic propagation of the rotation matrix, $\mathbf{A}$, as an
792   alternative method for the integration of orientational degrees of
793   freedom. The method starts with a straightforward splitting of the
794   Liouville operator:
795 < \begin{equation}
796 < eq here
797 < \label{introEq:SR1}
798 < \end{equation}
799 < Where $\boldsymbol{\tau}(\mathbf{A})$ are the torques of the system
800 < due to the configuration, and $\boldsymbol{/pi}$ are the conjugate
795 > \begin{align}
796 > iL_{\text{pos}} &= \dot{q}\frac{\partial}{\partial q} +
797 >        \mathbf{\dot{A}}\frac{\partial}{\partial \mathbf{A}}
798 > \label{introEq:SR1a} \\%
799 > %
800 > iL_F &= F(q)\frac{\partial}{\partial p}
801 > \label{introEq:SR1b} \\%
802 > iL_{\tau} &= \tau(\mathbf{A})\frac{\partial}{\partial \pi}
803 > \label{introEq:SR1b} \\%
804 > \end{align}
805 > Where $\tau(\mathbf{A})$ is the torque of the system
806 > due to the configuration, and $\pi$ is the conjugate
807   angular momenta of the system. The propagator, $G(\Delta t)$, becomes
808   \begin{equation}
809 < eq here
809 > G(\Delta t) = e^{\frac{\Delta t}{2} F(q)\frac{\partial}{\partial p}} \,
810 >        e^{\frac{\Delta t}{2} \tau(\mathbf{A})\frac{\partial}{\partial \pi}} \,
811 >        e^{\Delta t\,iL_{\text{pos}}} \,
812 >        e^{\frac{\Delta t}{2} \tau(\mathbf{A})\frac{\partial}{\partial \pi}} \,
813 >        e^{\frac{\Delta t}{2} F(q)\frac{\partial}{\partial p}}
814   \label{introEq:SR2}
815   \end{equation}
816   Propagation of the linear and angular momenta follows as in the Verlet
817   scheme.  The propagation of positions also follows the Verlet scheme
818   with the addition of a further symplectic splitting of the rotation
819 < matrix propagation, $\mathcal{G}_{\text{rot}}(\Delta t)$.
819 > matrix propagation, $\mathcal{U}_{\text{rot}}(\Delta t)$, within
820 > $U_{\text{pos}}(\Delta t)$.
821   \begin{equation}
822 < eq here
822 > \mathcal{U}_{\text{rot}}(\Delta t) =
823 >        \mathcal{U}_x \biggl(\frac{\Delta t}{2}\biggr)\,
824 >        \mathcal{U}_y \biggl(\frac{\Delta t}{2}\biggr)\,
825 >        \mathcal{U}_z (\Delta t)\,
826 >        \mathcal{U}_y \biggl(\frac{\Delta t}{2}\biggr)\,
827 >        \mathcal{U}_x \biggl(\frac{\Delta t}{2}\biggr)\,
828   \label{introEq:SR3}
829   \end{equation}
830 < Where $\mathcal{G}_j$ is a unitary rotation of $\mathbf{A}$ and
831 < $\boldsymbol{\pi}$ about each axis $j$.  As all propagations are now
830 > Where $\mathcal{U}_j$ is a unitary rotation of $\mathbf{A}$ and
831 > $\pi$ about each axis $j$.  As all propagations are now
832   unitary and symplectic, the entire integration scheme is also
833   symplectic and time reversible.
834  
835   \section{\label{introSec:layout}Dissertation Layout}
836  
837 < This dissertation is divided as follows:Chapt.~\ref{chapt:RSA}
837 > This dissertation is divided as follows:Ch.~\ref{chapt:RSA}
838   presents the random sequential adsorption simulations of related
839   pthalocyanines on a gold (111) surface. Ch.~\ref{chapt:OOPSE}
840   is about the writing of the molecular dynamics simulation package
841 < {\sc oopse}, Ch.~\ref{chapt:lipid} regards the simulations of
842 < phospholipid bilayers using a mesoscale model, and lastly,
841 > {\sc oopse}. Ch.~\ref{chapt:lipid} regards the simulations of
842 > phospholipid bilayers using a mesoscale model. And lastly,
843   Ch.~\ref{chapt:conclusion} concludes this dissertation with a
844   summary of all results. The chapters are arranged in chronological
845   order, and reflect the progression of techniques I employed during my
846   research.  
847  
848 < The chapter concerning random sequential adsorption
849 < simulations is a study in applying the principles of theoretical
850 < research in order to obtain a simple model capable of explaining the
851 < results.  My advisor, Dr. Gezelter, and I were approached by a
852 < colleague, Dr. Lieberman, about possible explanations for partial
853 < coverage of a gold surface by a particular compound of hers. We
854 < suggested it might be due to the statistical packing fraction of disks
855 < on a plane, and set about to simulate this system.  As the events in
856 < our model were not dynamic in nature, a Monte Carlo method was
857 < employed.  Here, if a molecule landed on the surface without
858 < overlapping another, then its landing was accepted.  However, if there
859 < was overlap, the landing we rejected and a new random landing location
860 < was chosen.  This defined our acceptance rules and allowed us to
861 < construct a Markov chain whose limiting distribution was the surface
862 < coverage in which we were interested.
848 > The chapter concerning random sequential adsorption simulations is a
849 > study in applying Statistical Mechanics simulation techniques in order
850 > to obtain a simple model capable of explaining the results.  My
851 > advisor, Dr. Gezelter, and I were approached by a colleague,
852 > Dr. Lieberman, about possible explanations for the  partial coverage of a
853 > gold surface by a particular compound of hers. We suggested it might
854 > be due to the statistical packing fraction of disks on a plane, and
855 > set about to simulate this system.  As the events in our model were
856 > not dynamic in nature, a Monte Carlo method was employed.  Here, if a
857 > molecule landed on the surface without overlapping another, then its
858 > landing was accepted.  However, if there was overlap, the landing we
859 > rejected and a new random landing location was chosen.  This defined
860 > our acceptance rules and allowed us to construct a Markov chain whose
861 > limiting distribution was the surface coverage in which we were
862 > interested.
863  
864   The following chapter, about the simulation package {\sc oopse},
865   describes in detail the large body of scientific code that had to be
866 < written in order to study phospholipid bilayer.  Although there are
866 > written in order to study phospholipid bilayers.  Although there are
867   pre-existing molecular dynamic simulation packages available, none
868   were capable of implementing the models we were developing.{\sc oopse}
869   is a unique package capable of not only integrating the equations of
# Line 809 | Line 875 | This model retains information about solvent ordering
875  
876   Bringing us to Ch.~\ref{chapt:lipid}. Using {\sc oopse}, I have been
877   able to parameterize a mesoscale model for phospholipid simulations.
878 < This model retains information about solvent ordering about the
878 > This model retains information about solvent ordering around the
879   bilayer, as well as information regarding the interaction of the
880 < phospholipid head groups' dipole with each other and the surrounding
880 > phospholipid head groups' dipoles with each other and the surrounding
881   solvent.  These simulations give us insight into the dynamic events
882   that lead to the formation of phospholipid bilayers, as well as
883   provide the foundation for future exploration of bilayer phase

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines