ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/xDissertation/Introduction.tex
Revision: 3354
Committed: Sat Mar 1 22:11:41 2008 UTC (16 years, 6 months ago) by xsun
Content type: application/x-tex
File size: 29712 byte(s)
Log Message:
writing up the dissertation.

File Contents

# User Rev Content
1 xsun 3336 \chapter{\label{chap:intro}INTRODUCTION AND BACKGROUND}
2 xsun 3347
3     \section{Background on the Problem\label{In:sec:pro}}
4     Phospholipid molecules are chosen to be studied in this dissertation
5     because of their critical role as a foundation of the bio-membrane
6     construction. The self assembled bilayer of the lipids when dispersed
7     in water is the micro structure of the membrane. The phase behavior of
8     lipid bilayer is explored experimentally~\cite{Cevc87}, however, fully
9     understanding on the mechanism is far beyond accomplished.
10    
11     \subsection{Ripple Phase\label{In:ssec:ripple}}
12     The {\it ripple phase} $P_{\beta'}$ of lipid bilayers, named from the
13     periodic buckling of the membrane, is an intermediate phase which is
14     developed either from heating the gel phase $L_{\beta'}$ or cooling
15     the fluid phase $L_\alpha$. Although the ripple phase is observed in
16     different
17     experiments~\cite{Sun96,Katsaras00,Copeland80,Meyer96,Kaasgaard03},
18     the mechanism of the formation of the ripple phase has never been
19     explained and the microscopic structure of the ripple phase has never
20     been elucidated by experiments. Computational simulation is a perfect
21     tool to study the microscopic properties for a system, however, the
22     long range dimension of the ripple structure and the long time scale
23     of the formation of the ripples are crucial obstacles to performing
24     the actual work. The idea to break through this dilemma forks into:
25     \begin{itemize}
26     \item Simplify the lipid model.
27     \item Improve the integrating algorithm.
28     \end{itemize}
29     In Ch.~\ref{chap:mc} and~\ref{chap:md}, we use a simple point dipole
30     model and a coarse-grained model to perform the Monte Carlo and
31     Molecular Dynamics simulations respectively, and in Ch.~\ref{chap:ld},
32     we implement a Langevin Dynamics algorithm to exclude the explicit
33     solvent to improve the efficiency of the simulations.
34    
35     \subsection{Lattice Model\label{In:ssec:model}}
36     The gel like characteristic of the ripple phase ensures the feasiblity
37     of applying the lattice model to study the system. It is claimed that
38     the packing of the lipid molecules in ripple phase is
39     hexagonal~\cite{Cevc87}. The popular $2$ dimensional lattice models,
40     {\it i.e.}, Ising model, Heisenberg model and $X-Y$ model, show
41     {\it frustration} on triangular lattice.
42     \begin{figure}
43     \centering
44 xsun 3354 \includegraphics[width=\linewidth]{./figures/inFrustration.pdf}
45 xsun 3347 \caption{Sketch to illustrate the frustration on triangular
46     lattice. Spins are represented by arrows, no matter which direction
47     the spin on the top of triangle points to, the Hamiltonian of the
48     system is the same, hence there are infinite possibilities for the
49     packing of the spins.}
50 xsun 3354 \label{Infig:frustration}
51 xsun 3347 \end{figure}
52 xsun 3354 Figure~\ref{Infig:frustration} shows an illustration of the frustration
53 xsun 3347 on a triangular lattice. The direction of the spin on top of the
54     triangle has no effects on the Hamiltonian of the system, therefore
55     infinite possibilities for the packing of spins induce the frustration
56     of the lattice.
57    
58     The lack of translational degree of freedom in lattice models prevents
59     their utilization on investigating the emergence of the surface
60     buckling which is the imposition of the ripple formation. In this
61     dissertation, a modified lattice model is introduced to this specific
62     situation in Ch.~\ref{chap:mc}.
63    
64     \section{Overview of Classical Statistical Mechanics\label{In:sec:SM}}
65     Statistical mechanics provides a way to calculate the macroscopic
66     properties of a system from the molecular interactions used in
67     computational simulations. This section serves as a brief introduction
68     to key concepts of classical statistical mechanics that we used in
69     this dissertation. Tolman gives an excellent introduction to the
70     principles of statistical mechanics~\cite{Tolman1979}. A large part of
71     section~\ref{In:sec:SM} will follow Tolman's notation.
72    
73     \subsection{Ensembles\label{In:ssec:ensemble}}
74     In classical mechanics, the state of the system is completely
75     described by the positions and momenta of all particles. If we have an
76     $N$ particle system, there are $6N$ coordinates ($3N$ position $(q_1,
77     q_2, \ldots, q_{3N})$ and $3N$ momenta $(p_1, p_2, \ldots, p_{3N})$)
78     to define the instantaneous state of the system. Each single set of
79     the $6N$ coordinates can be considered as a unique point in a $6N$
80     dimensional space where each perpendicular axis is one of
81     $q_{i\alpha}$ or $p_{i\alpha}$ ($i$ is the particle and $\alpha$ is
82     the spatial axis). This $6N$ dimensional space is known as phase
83     space. The instantaneous state of the system is a single point in
84     phase space. A trajectory is a connected path of points in phase
85     space. An ensemble is a collection of systems described by the same
86     macroscopic observables but which have microscopic state
87     distributions. In phase space an ensemble is a collection of a set of
88     representive points. A density distribution $\rho(q^N, p^N)$ of the
89     representive points in phase space describes the condition of an
90     ensemble of identical systems. Since this density may change with
91     time, it is also a function of time. $\rho(q^N, p^N, t)$ describes the
92     ensemble at a time $t$, and $\rho(q^N, p^N, t')$ describes the
93     ensemble at a later time $t'$. For convenience, we will use $\rho$
94     instead of $\rho(q^N, p^N, t)$ in the following disccusion. If we
95     normalize $\rho$ to unity,
96     \begin{equation}
97     1 = \int d \vec q~^N \int d \vec p~^N \rho,
98 xsun 3354 \label{Ineq:normalized}
99 xsun 3347 \end{equation}
100     then the value of $\rho$ gives the probability of finding the system
101     in a unit volume in the phase space.
102    
103     Liouville's theorem describes the change in density $\rho$ with
104     time. The number of representive points at a given volume in the phase
105     space at any instant can be written as:
106     \begin{equation}
107 xsun 3354 \label{Ineq:deltaN}
108 xsun 3347 \delta N = \rho~\delta q_1 \delta q_2 \ldots \delta q_N \delta p_1 \delta p_2 \ldots \delta p_N.
109     \end{equation}
110     To calculate the change in the number of representive points in this
111     volume, let us consider a simple condition: the change in the number
112     of representive points in $q_1$ axis. The rate of the number of the
113     representive points entering the volume at $q_1$ per unit time is:
114     \begin{equation}
115 xsun 3354 \label{Ineq:deltaNatq1}
116 xsun 3347 \rho~\dot q_1 \delta q_2 \ldots \delta q_N \delta p_1 \delta p_2 \ldots \delta p_N,
117     \end{equation}
118     and the rate of the number of representive points leaving the volume
119     at another position $q_1 + \delta q_1$ is:
120     \begin{equation}
121 xsun 3354 \label{Ineq:deltaNatq1plusdeltaq1}
122 xsun 3347 \left( \rho + \frac{\partial \rho}{\partial q_1} \delta q_1 \right)\left(\dot q_1 +
123     \frac{\partial \dot q_1}{\partial q_1} \delta q_1 \right)\delta q_2 \ldots \delta q_N \delta p_1 \delta p_2 \ldots \delta p_N.
124     \end{equation}
125     Here the higher order differentials are neglected. So the change of
126     the number of the representive points is the difference of
127 xsun 3354 eq.~\ref{Ineq:deltaNatq1} and eq.~\ref{Ineq:deltaNatq1plusdeltaq1}, which
128 xsun 3347 gives us:
129     \begin{equation}
130 xsun 3354 \label{Ineq:deltaNatq1axis}
131 xsun 3347 -\left(\rho \frac{\partial {\dot q_1}}{\partial q_1} + \frac{\partial {\rho}}{\partial q_1} \dot q_1 \right)\delta q_1 \delta q_2 \ldots \delta q_N \delta p_1 \delta p_2 \ldots \delta p_N,
132     \end{equation}
133     where, higher order differetials are neglected. If we sum over all the
134     axes in the phase space, we can get the change of the number of
135     representive points in a given volume with time:
136     \begin{equation}
137 xsun 3354 \label{Ineq:deltaNatGivenVolume}
138 xsun 3347 \frac{d(\delta N)}{dt} = -\sum_{i=1}^N \left[\rho \left(\frac{\partial
139     {\dot q_i}}{\partial q_i} + \frac{\partial
140     {\dot p_i}}{\partial p_i}\right) + \left( \frac{\partial {\rho}}{\partial
141     q_i} \dot q_i + \frac{\partial {\rho}}{\partial p_i} \dot p_i\right)\right]\delta q_1 \delta q_2 \ldots \delta q_N \delta p_1 \delta p_2 \ldots \delta p_N.
142     \end{equation}
143     From Hamilton's equation of motion,
144     \begin{equation}
145     \frac{\partial {\dot q_i}}{\partial q_i} = - \frac{\partial
146     {\dot p_i}}{\partial p_i},
147 xsun 3354 \label{Ineq:canonicalFormOfEquationOfMotion}
148 xsun 3347 \end{equation}
149     this cancels out the first term on the right side of
150 xsun 3354 eq.~\ref{Ineq:deltaNatGivenVolume}. If both sides of
151     eq.~\ref{Ineq:deltaNatGivenVolume} are divided by $\delta q_1 \delta q_2
152 xsun 3347 \ldots \delta q_N \delta p_1 \delta p_2 \ldots \delta p_N$, then we
153     can derive Liouville's theorem:
154     \begin{equation}
155     \left( \frac{\partial \rho}{\partial t} \right)_{q, p} = -\sum_{i} \left(
156     \frac{\partial {\rho}}{\partial
157     q_i} \dot q_i + \frac{\partial {\rho}}{\partial p_i} \dot p_i \right).
158 xsun 3354 \label{Ineq:simpleFormofLiouville}
159 xsun 3347 \end{equation}
160     This is the basis of statistical mechanics. If we move the right
161 xsun 3354 side of equation~\ref{Ineq:simpleFormofLiouville} to the left, we
162 xsun 3347 will obtain
163     \begin{equation}
164     \left( \frac{\partial \rho}{\partial t} \right)_{q, p} + \sum_{i} \left(
165     \frac{\partial {\rho}}{\partial
166     q_i} \dot q_i + \frac{\partial {\rho}}{\partial p_i} \dot p_i \right)
167     = 0.
168 xsun 3354 \label{Ineq:anotherFormofLiouville}
169 xsun 3347 \end{equation}
170     It is easy to note that the left side of
171 xsun 3354 equation~\ref{Ineq:anotherFormofLiouville} is the total derivative of
172 xsun 3347 $\rho$ with respect of $t$, which means
173     \begin{equation}
174     \frac{d \rho}{dt} = 0,
175 xsun 3354 \label{Ineq:conservationofRho}
176 xsun 3347 \end{equation}
177     and the rate of density change is zero in the neighborhood of any
178     selected moving representive points in the phase space.
179    
180     The condition of the ensemble is determined by the density
181     distribution. If we consider the density distribution as only a
182     function of $q$ and $p$, which means the rate of change of the phase
183     space density in the neighborhood of all representive points in the
184     phase space is zero,
185     \begin{equation}
186     \left( \frac{\partial \rho}{\partial t} \right)_{q, p} = 0.
187 xsun 3354 \label{Ineq:statEquilibrium}
188 xsun 3347 \end{equation}
189     We may conclude the ensemble is in {\it statistical equilibrium}. An
190     ensemble in statistical equilibrium often means the system is also in
191     macroscopic equilibrium. If $\left( \frac{\partial \rho}{\partial t}
192     \right)_{q, p} = 0$, then
193     \begin{equation}
194     \sum_{i} \left(
195     \frac{\partial {\rho}}{\partial
196     q_i} \dot q_i + \frac{\partial {\rho}}{\partial p_i} \dot p_i \right)
197     = 0.
198 xsun 3354 \label{Ineq:constantofMotion}
199 xsun 3347 \end{equation}
200     If $\rho$ is a function only of some constant of the motion, $\rho$ is
201     independent of time. For a conservative system, the energy of the
202     system is one of the constants of the motion. Here are several
203     examples: when the density distribution is constant everywhere in the
204     phase space,
205     \begin{equation}
206     \rho = \mathrm{const.}
207 xsun 3354 \label{Ineq:uniformEnsemble}
208 xsun 3347 \end{equation}
209     the ensemble is called {\it uniform ensemble}. Another useful
210     ensemble is called {\it microcanonical ensemble}, for which:
211     \begin{equation}
212     \rho = \delta \left( H(q^N, p^N) - E \right) \frac{1}{\Sigma (N, V, E)}
213 xsun 3354 \label{Ineq:microcanonicalEnsemble}
214 xsun 3347 \end{equation}
215     where $\Sigma(N, V, E)$ is a normalization constant parameterized by
216     $N$, the total number of particles, $V$, the total physical volume and
217     $E$, the total energy. The physical meaning of $\Sigma(N, V, E)$ is
218     the phase space volume accessible to a microcanonical system with
219     energy $E$ evolving under Hamilton's equations. $H(q^N, p^N)$ is the
220     Hamiltonian of the system. The Gibbs entropy is defined as
221     \begin{equation}
222     S = - k_B \int d \vec q~^N \int d \vec p~^N \rho \ln [C^N \rho],
223 xsun 3354 \label{Ineq:gibbsEntropy}
224 xsun 3347 \end{equation}
225     where $k_B$ is the Boltzmann constant and $C^N$ is a number which
226     makes the argument of $\ln$ dimensionless, in this case, it is the
227     total phase space volume of one state. The entropy in microcanonical
228     ensemble is given by
229     \begin{equation}
230     S = k_B \ln \left(\frac{\Sigma(N, V, E)}{C^N}\right).
231 xsun 3354 \label{Ineq:entropy}
232 xsun 3347 \end{equation}
233     If the density distribution $\rho$ is given by
234     \begin{equation}
235     \rho = \frac{1}{Z_N}e^{-H(q^N, p^N) / k_B T},
236 xsun 3354 \label{Ineq:canonicalEnsemble}
237 xsun 3347 \end{equation}
238     the ensemble is known as the {\it canonical ensemble}. Here,
239     \begin{equation}
240     Z_N = \int d \vec q~^N \int_\Gamma d \vec p~^N e^{-H(q^N, p^N) / k_B T},
241 xsun 3354 \label{Ineq:partitionFunction}
242 xsun 3347 \end{equation}
243     which is also known as {\it partition function}. $\Gamma$ indicates
244     that the integral is over all the phase space. In the canonical
245     ensemble, $N$, the total number of particles, $V$, total volume, and
246     $T$, the temperature are constants. The systems with the lowest
247     energies hold the largest population. According to maximum principle,
248     the thermodynamics maximizes the entropy $S$,
249     \begin{equation}
250     \begin{array}{ccc}
251     \delta S = 0 & \mathrm{and} & \delta^2 S < 0.
252     \end{array}
253 xsun 3354 \label{Ineq:maximumPrinciple}
254 xsun 3347 \end{equation}
255 xsun 3354 From Eq.~\ref{Ineq:maximumPrinciple} and two constrains of the canonical
256 xsun 3347 ensemble, {\it i.e.}, total probability and average energy conserved,
257     the partition function is calculated as
258     \begin{equation}
259     Z_N = e^{-A/k_B T},
260 xsun 3354 \label{Ineq:partitionFunctionWithFreeEnergy}
261 xsun 3347 \end{equation}
262     where $A$ is the Helmholtz free energy. The significance of
263 xsun 3354 Eq.~\ref{Ineq:entropy} and~\ref{Ineq:partitionFunctionWithFreeEnergy} is
264 xsun 3347 that they serve as a connection between macroscopic properties of the
265     system and the distribution of the microscopic states.
266    
267     There is an implicit assumption that our arguments are based on so
268     far. A representive point in the phase space is equally to be found in
269     any same extent of the regions. In other words, all energetically
270     accessible states are represented equally, the probabilities to find
271     the system in any of the accessible states is equal. This is called
272     equal a {\it priori} probabilities.
273    
274     \subsection{Statistical Average\label{In:ssec:average}}
275     Given the density distribution $\rho$ in the phase space, the average
276     of any quantity ($F(q^N, p^N$)) which depends on the coordinates
277     ($q^N$) and the momenta ($p^N$) for all the systems in the ensemble
278     can be calculated based on the definition shown by
279 xsun 3354 Eq.~\ref{Ineq:statAverage1}
280 xsun 3347 \begin{equation}
281     \langle F(q^N, p^N, t) \rangle = \frac{\int d \vec q~^N \int d \vec p~^N
282     F(q^N, p^N, t) \rho}{\int d \vec q~^N \int d \vec p~^N \rho}.
283 xsun 3354 \label{Ineq:statAverage1}
284 xsun 3347 \end{equation}
285     Since the density distribution $\rho$ is normalized to unity, the mean
286     value of $F(q^N, p^N)$ is simplified to
287     \begin{equation}
288     \langle F(q^N, p^N, t) \rangle = \int d \vec q~^N \int d \vec p~^N F(q^N,
289     p^N, t) \rho,
290 xsun 3354 \label{Ineq:statAverage2}
291 xsun 3347 \end{equation}
292     called {\it ensemble average}. However, the quantity is often averaged
293     for a finite time in real experiments,
294     \begin{equation}
295     \langle F(q^N, p^N) \rangle_t = \lim_{T \rightarrow \infty}
296     \frac{1}{T} \int_{t_0}^{t_0+T} F(q^N, p^N, t) dt.
297 xsun 3354 \label{Ineq:timeAverage1}
298 xsun 3347 \end{equation}
299     Usually this time average is independent of $t_0$ in statistical
300 xsun 3354 mechanics, so Eq.~\ref{Ineq:timeAverage1} becomes
301 xsun 3347 \begin{equation}
302     \langle F(q^N, p^N) \rangle_t = \lim_{T \rightarrow \infty}
303     \frac{1}{T} \int_{0}^{T} F(q^N, p^N, t) dt
304 xsun 3354 \label{Ineq:timeAverage2}
305 xsun 3347 \end{equation}
306     for an infinite time interval.
307    
308     {\it ergodic hypothesis}, an important hypothesis from the statistical
309     mechanics point of view, states that the system will eventually pass
310     arbitrarily close to any point that is energetically accessible in
311     phase space. Mathematically, this leads to
312     \begin{equation}
313     \langle F(q^N, p^N, t) \rangle = \langle F(q^N, p^N) \rangle_t.
314 xsun 3354 \label{Ineq:ergodicity}
315 xsun 3347 \end{equation}
316 xsun 3354 Eq.~\ref{Ineq:ergodicity} validates the Monte Carlo method which we will
317 xsun 3347 discuss in section~\ref{In:ssec:mc}. An ensemble average of a quantity
318     can be related to the time average measured in the experiments.
319    
320     \subsection{Correlation Function\label{In:ssec:corr}}
321     Thermodynamic properties can be computed by equillibrium statistical
322     mechanics. On the other hand, {\it Time correlation function} is a
323     powerful method to understand the evolution of a dynamic system in
324     non-equillibrium statistical mechanics. Imagine a property $A(q^N,
325     p^N, t)$ as a function of coordinates $q^N$ and momenta $p^N$ has an
326     intial value at $t_0$, at a later time $t_0 + \tau$ this value is
327     changed. If $\tau$ is very small, the change of the value is minor,
328     and the later value of $A(q^N, p^N, t_0 +
329     \tau)$ is correlated to its initial value. Howere, when $\tau$ is large,
330     this correlation is lost. The correlation function is a measurement of
331     this relationship and is defined by~\cite{Berne90}
332     \begin{equation}
333     C(t) = \langle A(0)A(\tau) \rangle = \lim_{T \rightarrow \infty}
334     \frac{1}{T} \int_{0}^{T} dt A(t) A(t + \tau).
335 xsun 3354 \label{Ineq:autocorrelationFunction}
336 xsun 3347 \end{equation}
337 xsun 3354 Eq.~\ref{Ineq:autocorrelationFunction} is the correlation function of a
338 xsun 3347 single variable, called {\it autocorrelation function}. The defination
339     of the correlation function for two different variables is similar to
340     that of autocorrelation function, which is
341     \begin{equation}
342     C(t) = \langle A(0)B(\tau) \rangle = \lim_{T \rightarrow \infty}
343     \frac{1}{T} \int_{0}^{T} dt A(t) B(t + \tau),
344 xsun 3354 \label{Ineq:crosscorrelationFunction}
345 xsun 3347 \end{equation}
346     and called {\it cross correlation function}.
347    
348 xsun 3354 In section~\ref{In:ssec:average} we know from Eq.~\ref{Ineq:ergodicity}
349 xsun 3347 the relationship between time average and ensemble average. We can put
350     the correlation function in a classical mechanics form,
351     \begin{equation}
352     C(t) = \langle A(0)A(\tau) \rangle = \int d \vec q~^N \int d \vec p~^N A(t) A(t + \tau) \rho(q, p)
353 xsun 3354 \label{Ineq:autocorrelationFunctionCM}
354 xsun 3347 \end{equation}
355     and
356     \begin{equation}
357     C(t) = \langle A(0)B(\tau) \rangle = \int d \vec q~^N \int d \vec p~^N A(t) B(t + \tau)
358     \rho(q, p)
359 xsun 3354 \label{Ineq:crosscorrelationFunctionCM}
360 xsun 3347 \end{equation}
361     as autocorrelation function and cross correlation function
362     respectively. $\rho(q, p)$ is the density distribution at equillibrium
363     in phase space. In many cases, the correlation function decay is a
364     single exponential
365     \begin{equation}
366     C(t) \sim e^{-t / \tau_r},
367 xsun 3354 \label{Ineq:relaxation}
368 xsun 3347 \end{equation}
369     where $\tau_r$ is known as relaxation time which discribes the rate of
370     the decay.
371    
372     \section{Methodolody\label{In:sec:method}}
373     The simulations performed in this dissertation are branched into two
374     main catalog, Monte Carlo and Molecular Dynamics. There are two main
375     difference between Monte Carlo and Molecular Dynamics simulations. One
376     is that the Monte Carlo simulation is time independent, and Molecular
377     Dynamics simulation is time involved. Another dissimilar is that the
378     Monte Carlo is a stochastic process, the configuration of the system
379     is not determinated by its past, however, using Moleuclar Dynamics,
380     the system is propagated by Newton's equation of motion, the
381     trajectory of the system evolved in the phase space is determined. A
382     brief introduction of the two algorithms are given in
383     section~\ref{In:ssec:mc} and~\ref{In:ssec:md}. An extension of the
384     Molecular Dynamics, Langevin Dynamics, is introduced by
385     section~\ref{In:ssec:ld}.
386    
387     \subsection{Monte Carlo\label{In:ssec:mc}}
388     Monte Carlo algorithm was first introduced by Metropolis {\it et
389     al.}.~\cite{Metropolis53} Basic Monte Carlo algorithm is usually
390     applied to the canonical ensemble, a Boltzmann-weighted ensemble, in
391     which the $N$, the total number of particles, $V$, total volume, $T$,
392     temperature are constants. The average energy is given by substituding
393 xsun 3354 Eq.~\ref{Ineq:canonicalEnsemble} into Eq.~\ref{Ineq:statAverage2},
394 xsun 3347 \begin{equation}
395     \langle E \rangle = \frac{1}{Z_N} \int d \vec q~^N \int d \vec p~^N E e^{-H(q^N, p^N) / k_B T}.
396 xsun 3354 \label{Ineq:energyofCanonicalEnsemble}
397 xsun 3347 \end{equation}
398     So are the other properties of the system. The Hamiltonian is the
399     summation of Kinetic energy $K(p^N)$ as a function of momenta and
400     Potential energy $U(q^N)$ as a function of positions,
401     \begin{equation}
402     H(q^N, p^N) = K(p^N) + U(q^N).
403 xsun 3354 \label{Ineq:hamiltonian}
404 xsun 3347 \end{equation}
405     If the property $A$ is only a function of position ($ A = A(q^N)$),
406     the mean value of $A$ is reduced to
407     \begin{equation}
408     \langle A \rangle = \frac{\int d \vec q~^N \int d \vec p~^N A e^{-U(q^N) / k_B T}}{\int d \vec q~^N \int d \vec p~^N e^{-U(q^N) / k_B T}},
409 xsun 3354 \label{Ineq:configurationIntegral}
410 xsun 3347 \end{equation}
411     The kinetic energy $K(p^N)$ is factored out in
412 xsun 3354 Eq.~\ref{Ineq:configurationIntegral}. $\langle A
413 xsun 3347 \rangle$ is a configuration integral now, and the
414 xsun 3354 Eq.~\ref{Ineq:configurationIntegral} is equivalent to
415 xsun 3347 \begin{equation}
416     \langle A \rangle = \int d \vec q~^N A \rho(q^N).
417 xsun 3354 \label{Ineq:configurationAve}
418 xsun 3347 \end{equation}
419    
420     In a Monte Carlo simulation of canonical ensemble, the probability of
421     the system being in a state $s$ is $\rho_s$, the change of this
422     probability with time is given by
423     \begin{equation}
424     \frac{d \rho_s}{dt} = \sum_{s'} [ -w_{ss'}\rho_s + w_{s's}\rho_{s'} ],
425 xsun 3354 \label{Ineq:timeChangeofProb}
426 xsun 3347 \end{equation}
427     where $w_{ss'}$ is the tansition probability of going from state $s$
428     to state $s'$. Since $\rho_s$ is independent of time at equilibrium,
429     \begin{equation}
430     \frac{d \rho_{s}^{equilibrium}}{dt} = 0,
431 xsun 3354 \label{Ineq:equiProb}
432 xsun 3347 \end{equation}
433     which means $\sum_{s'} [ -w_{ss'}\rho_s + w_{s's}\rho_{s'} ]$ is $0$
434     for all $s'$. So
435     \begin{equation}
436     \frac{\rho_s^{equilibrium}}{\rho_{s'}^{equilibrium}} = \frac{w_{s's}}{w_{ss'}}.
437 xsun 3354 \label{Ineq:relationshipofRhoandW}
438 xsun 3347 \end{equation}
439     If
440     \begin{equation}
441     \frac{w_{s's}}{w_{ss'}} = e^{-(U_s - U_{s'}) / k_B T},
442 xsun 3354 \label{Ineq:conditionforBoltzmannStatistics}
443 xsun 3347 \end{equation}
444     then
445     \begin{equation}
446     \frac{\rho_s^{equilibrium}}{\rho_{s'}^{equilibrium}} = e^{-(U_s - U_{s'}) / k_B T}.
447 xsun 3354 \label{Ineq:satisfyofBoltzmannStatistics}
448 xsun 3347 \end{equation}
449 xsun 3354 Eq.~\ref{Ineq:satisfyofBoltzmannStatistics} implies that
450 xsun 3347 $\rho^{equilibrium}$ satisfies Boltzmann statistics. An algorithm,
451     shows how Monte Carlo simulation generates a transition probability
452 xsun 3354 governed by \ref{Ineq:conditionforBoltzmannStatistics}, is schemed as
453 xsun 3347 \begin{enumerate}
454 xsun 3354 \item\label{Initm:oldEnergy} Choose an particle randomly, calculate the energy.
455     \item\label{Initm:newEnergy} Make a random displacement for particle,
456 xsun 3347 calculate the new energy.
457     \begin{itemize}
458 xsun 3354 \item Keep the new configuration and return to step~\ref{Initm:oldEnergy} if energy
459 xsun 3347 goes down.
460     \item Pick a random number between $[0,1]$ if energy goes up.
461     \begin{itemize}
462     \item Keep the new configuration and return to
463 xsun 3354 step~\ref{Initm:oldEnergy} if the random number smaller than
464 xsun 3347 $e^{-(U_{new} - U_{old})} / k_B T$.
465     \item Keep the old configuration and return to
466 xsun 3354 step~\ref{Initm:oldEnergy} if the random number larger than
467 xsun 3347 $e^{-(U_{new} - U_{old})} / k_B T$.
468     \end{itemize}
469     \end{itemize}
470 xsun 3354 \item\label{Initm:accumulateAvg} Accumulate the average after it converges.
471 xsun 3347 \end{enumerate}
472     It is important to notice that the old configuration has to be sampled
473     again if it is kept.
474    
475     \subsection{Molecular Dynamics\label{In:ssec:md}}
476     Although some of properites of the system can be calculated from the
477     ensemble average in Monte Carlo simulations, due to the nature of
478     lacking in the time dependence, it is impossible to gain information
479     of those dynamic properties from Monte Carlo simulations. Molecular
480     Dynamics is a measurement of the evolution of the positions and
481     momenta of the particles in the system. The evolution of the system
482     obeys laws of classical mechanics, in most situations, there is no
483     need for the count of the quantum effects. For a real experiment, the
484     instantaneous positions and momenta of the particles in the system are
485     neither important nor measurable, the observable quantities are
486     usually a average value for a finite time interval. These quantities
487     are expressed as a function of positions and momenta in Melecular
488     Dynamics simulations. Like the thermal temperature of the system is
489     defined by
490     \begin{equation}
491     \frac{1}{2} k_B T = \langle \frac{1}{2} m v_\alpha \rangle,
492 xsun 3354 \label{Ineq:temperature}
493 xsun 3347 \end{equation}
494     here $m$ is the mass of the particle and $v_\alpha$ is the $\alpha$
495     component of the velocity of the particle. The right side of
496 xsun 3354 Eq.~\ref{Ineq:temperature} is the average kinetic energy of the
497 xsun 3347 system. A simple Molecular Dynamics simulation scheme
498     is:~\cite{Frenkel1996}
499     \begin{enumerate}
500 xsun 3354 \item\label{Initm:initialize} Assign the initial positions and momenta
501 xsun 3347 for the particles in the system.
502 xsun 3354 \item\label{Initm:calcForce} Calculate the forces.
503     \item\label{Initm:equationofMotion} Integrate the equation of motion.
504 xsun 3347 \begin{itemize}
505 xsun 3354 \item Return to step~\ref{Initm:calcForce} if the equillibrium is
506 xsun 3347 not achieved.
507 xsun 3354 \item Go to step~\ref{Initm:calcAvg} if the equillibrium is
508 xsun 3347 achieved.
509     \end{itemize}
510 xsun 3354 \item\label{Initm:calcAvg} Compute the quantities we are interested in.
511 xsun 3347 \end{enumerate}
512     The initial positions of the particles are chosen as that there is no
513     overlap for the particles. The initial velocities at first are
514     distributed randomly to the particles, and then shifted to make the
515     momentum of the system $0$, at last scaled to the desired temperature
516 xsun 3354 of the simulation according Eq.~\ref{Ineq:temperature}.
517 xsun 3347
518 xsun 3354 The core of Molecular Dynamics simulations is step~\ref{Initm:calcForce}
519     and~\ref{Initm:equationofMotion}. The calculation of the forces are
520 xsun 3347 often involved numerous effort, this is the most time consuming step
521     in the Molecular Dynamics scheme. The evaluation of the forces is
522     followed by
523     \begin{equation}
524     f(q) = - \frac{\partial U(q)}{\partial q},
525 xsun 3354 \label{Ineq:force}
526 xsun 3347 \end{equation}
527     $U(q)$ is the potential of the system. Once the forces computed, are
528     the positions and velocities updated by integrating Newton's equation
529     of motion,
530     \begin{equation}
531     f(q) = \frac{dp}{dt} = \frac{m dv}{dt}.
532 xsun 3354 \label{Ineq:newton}
533 xsun 3347 \end{equation}
534     Here is an example of integrating algorithms, Verlet algorithm, which
535     is one of the best algorithms to integrate Newton's equation of
536     motion. The Taylor expension of the position at time $t$ is
537     \begin{equation}
538     q(t+\Delta t)= q(t) + v(t) \Delta t + \frac{f(t)}{2m}\Delta t^2 +
539     \frac{\Delta t^3}{3!}\frac{\partial^3 q(t)}{\partial t^3} +
540     \mathcal{O}(\Delta t^4)
541 xsun 3354 \label{Ineq:verletFuture}
542 xsun 3347 \end{equation}
543     for a later time $t+\Delta t$, and
544     \begin{equation}
545     q(t-\Delta t)= q(t) - v(t) \Delta t + \frac{f(t)}{2m}\Delta t^2 -
546     \frac{\Delta t^3}{3!}\frac{\partial^3 q(t)}{\partial t^3} +
547     \mathcal{O}(\Delta t^4) ,
548 xsun 3354 \label{Ineq:verletPrevious}
549 xsun 3347 \end{equation}
550     for a previous time $t-\Delta t$. The summation of the
551 xsun 3354 Eq.~\ref{Ineq:verletFuture} and~\ref{Ineq:verletPrevious} gives
552 xsun 3347 \begin{equation}
553     q(t+\Delta t)+q(t-\Delta t) =
554     2q(t) + \frac{f(t)}{m}\Delta t^2 + \mathcal{O}(\Delta t^4),
555 xsun 3354 \label{Ineq:verletSum}
556 xsun 3347 \end{equation}
557     so, the new position can be expressed as
558     \begin{equation}
559     q(t+\Delta t) \approx
560     2q(t) - q(t-\Delta t) + \frac{f(t)}{m}\Delta t^2.
561 xsun 3354 \label{Ineq:newPosition}
562 xsun 3347 \end{equation}
563     The higher order of the $\Delta t$ is omitted.
564    
565     Numerous technics and tricks are applied to Molecular Dynamics
566     simulation to gain more efficiency or more accuracy. The simulation
567     engine used in this dissertation for the Molecular Dynamics
568     simulations is {\sc oopse}, more details of the algorithms and
569     technics can be found in~\cite{Meineke2005}.
570    
571     \subsection{Langevin Dynamics\label{In:ssec:ld}}
572     In many cases, the properites of the solvent in a system, like the
573     lipid-water system studied in this dissertation, are less important to
574     the researchers. However, the major computational expense is spent on
575     the solvent in the Molecular Dynamics simulations because of the large
576     number of the solvent molecules compared to that of solute molecules,
577     the ratio of the number of lipid molecules to the number of water
578     molecules is $1:25$ in our lipid-water system. The efficiency of the
579     Molecular Dynamics simulations is greatly reduced.
580    
581     As an extension of the Molecular Dynamics simulations, the Langevin
582     Dynamics seeks a way to avoid integrating equation of motion for
583     solvent particles without losing the Brownian properites of solute
584     particles. A common approximation is that the coupling of the solute
585     and solvent is expressed as a set of harmonic oscillators. So the
586     Hamiltonian of such a system is written as
587     \begin{equation}
588     H = \frac{p^2}{2m} + U(q) + H_B + \Delta U(q),
589 xsun 3354 \label{Ineq:hamiltonianofCoupling}
590 xsun 3347 \end{equation}
591     where $H_B$ is the Hamiltonian of the bath which equals to
592     \begin{equation}
593     H_B = \sum_{\alpha = 1}^{N} \left\{ \frac{p_\alpha^2}{2m_\alpha} +
594     \frac{1}{2} m_\alpha \omega_\alpha^2 q_\alpha^2\right\},
595 xsun 3354 \label{Ineq:hamiltonianofBath}
596 xsun 3347 \end{equation}
597     $\alpha$ is all the degree of freedoms of the bath, $\omega$ is the
598     bath frequency, and $\Delta U(q)$ is the bilinear coupling given by
599     \begin{equation}
600     \Delta U = -\sum_{\alpha = 1}^{N} g_\alpha q_\alpha q,
601 xsun 3354 \label{Ineq:systemBathCoupling}
602 xsun 3347 \end{equation}
603     where $g$ is the coupling constant. By solving the Hamilton's equation
604     of motion, the {\it Generalized Langevin Equation} for this system is
605     derived to
606     \begin{equation}
607     m \ddot q = -\frac{\partial W(q)}{\partial q} - \int_0^t \xi(t) \dot q(t-t')dt' + R(t),
608 xsun 3354 \label{Ineq:gle}
609 xsun 3347 \end{equation}
610     with mean force,
611     \begin{equation}
612     W(q) = U(q) - \sum_{\alpha = 1}^N \frac{g_\alpha^2}{2m_\alpha
613     \omega_\alpha^2}q^2,
614 xsun 3354 \label{Ineq:meanForce}
615 xsun 3347 \end{equation}
616     being only a dependence of coordinates of the solute particles, {\it
617     friction kernel},
618     \begin{equation}
619     \xi(t) = \sum_{\alpha = 1}^N \frac{-g_\alpha}{m_\alpha
620     \omega_\alpha} \cos(\omega_\alpha t),
621 xsun 3354 \label{Ineq:xiforGLE}
622 xsun 3347 \end{equation}
623     and the random force,
624     \begin{equation}
625     R(t) = \sum_{\alpha = 1}^N \left( g_\alpha q_\alpha(0)-\frac{g_\alpha}{m_\alpha
626     \omega_\alpha^2}q(0)\right) \cos(\omega_\alpha t) + \frac{\dot
627     q_\alpha(0)}{\omega_\alpha} \sin(\omega_\alpha t),
628 xsun 3354 \label{Ineq:randomForceforGLE}
629 xsun 3347 \end{equation}
630     as only a dependence of the initial conditions. The relationship of
631     friction kernel $\xi(t)$ and random force $R(t)$ is given by
632     \begin{equation}
633     \xi(t) = \frac{1}{k_B T} \langle R(t)R(0) \rangle
634 xsun 3354 \label{Ineq:relationshipofXiandR}
635 xsun 3347 \end{equation}
636     from their definitions. In Langevin limit, the friction is treated
637     static, which means
638     \begin{equation}
639     \xi(t) = 2 \xi_0 \delta(t).
640 xsun 3354 \label{Ineq:xiofStaticFriction}
641 xsun 3347 \end{equation}
642 xsun 3354 After substitude $\xi(t)$ into Eq.~\ref{Ineq:gle} with
643     Eq.~\ref{Ineq:xiofStaticFriction}, {\it Langevin Equation} is extracted
644 xsun 3347 to
645     \begin{equation}
646     m \ddot q = -\frac{\partial U(q)}{\partial q} - \xi \dot q(t) + R(t).
647 xsun 3354 \label{Ineq:langevinEquation}
648 xsun 3347 \end{equation}
649     The applying of Langevin Equation to dynamic simulations is discussed
650     in Ch.~\ref{chap:ld}.