\chapter{Molecular Simulations}
\label{ch:methods}
Molecular simulation is a generic term encompassing both Monte Carlo (MC) and Molecular Dynamics (MD) computing methods. These methods can be used to provide exact results for problems in statistical mechanics which would otherwise only be solved by approximate methods, or are intractable. Molecular simulations can provide a direct route from microscopic details of a system to macroscopic properties of experimental interest.  It is also used to obtain information about processes that are difficult to probe experimentally.
\par
The MD simulation technique follows the natural time evolution of the system, which allows the prediction of static and dynamic properties of the system directly from  underlying interactions between molecules. However, sampling problems in MD are connected to the difficulty of spanning broad scales of time and this problem is hard to surmount because of the need to adhere to the governing dynamics.  On the other hand, the MC method uses random number generators to sample the phase space and  clever algorithms can be developed to optimize the efficiency of the sampling.
\section{Theory}
The microscopic state of a system can be specified in terms of the position and momentum of the constituent set of particles. If we make the approximation that a classical description is adequate, the Hamiltonian $H$ of a system of $N$ particles can be written as a sum of the kinetic and potential energy functions of the set of coordinates $\mathbf q_i$ and momenta $\mathbf p_i$ of each particle $i$:\cite{Allen1989}
%
\begin{equation}
H_{(\mathbf{q,p})}= K_{(\mathbf p)} + U_{(\mathbf q)},
\end{equation}
$\mathbf q$ and $\mathbf p$ are generalized coordinates and conjugate momenta, respectively. To obtain thermodynamic averages over an ensemble, which is characterized by the macroscopic variables $(N, V, T)$, where $N$ is the number of particles, $V$ is the volume and $T$ is the temperature of the system, it is necessary to know the probability of finding the system at each and every point in the phase space. This probability distribution is given by the Boltzmann distribution function
%
\begin{equation}
\label{bprob}
\rho (\mathbf{q,p}) = \frac{exp\left[{\frac{-H(\mathbf{q,p})}{k_BT}}\right]}{Z},
\end{equation}
%
where the canonical partition function $Z$ is the integral over the phase space:
%
\begin{equation}
Z = \int d\mathbf q\int d\mathbf p\ exp\left[{\frac{-H(\mathbf{q,p})}{k_BT}}\right].
\end{equation}
%
Once this distribution function is known, it can be used to calculate the phase space averages of any dynamic variable $A(\mathbf{q,p})$ of interest:
%
\begin{equation}
\langle A(\mathbf{q,p} )\rangle =\int d\mathbf q\int d\mathbf p\ \rho (\mathbf{q,p}) A(\mathbf{q,p}).
\end{equation}
%
However, in order to calculate these thermodynamic averages, it is necessary to simultaneously know the Boltzmann probability (Eq. \ref{bprob}) for each and every state $\{\mathbf{q,p}\}$, which is not feasible for many particle systems. An alternative strategy for calculating system wide averages is to follow the motion of a single point through the phase space instead of averaging over the whole phase space at once, i.e., the motion of a single point (a single molecular state) through the phase space is followed as a function of time, and the averages are calculated only over those points that were visited during the excursion. Averages calculated in this way are called ``dynamic averages''. The motion of a single point through the phase space is obtained by integrating the equation of motion of the system.\cite{Oren2001}
Dynamic averages of any dynamical variable $A(\mathbf{q,p} )$ can now be calculated along this trajectory as follows:
%
\begin{equation}
\langle A(\mathbf{q,p} )\rangle = \frac{1}{\tau}\int\limits^\tau _0 A\Big( \mathbf q (t), \mathbf p(t)\Big)dt,
\end{equation}
%
where $\tau$ is the duration of the simulation. It is hoped that the point that is being dynamically followed will eventually cover all the phase space and that the dynamic average will converge to the desired thermodynamic average. A key concept that ties the two averaging strategies together is the \emph{ergodic hypothesis}. This hypothesis states that, for an infinitely long trajectory, the thermodynamic ensemble average and the dynamic average become equivalent to each other.
%
\section{Force fields}
The core of any force field is the potential energy function used to describe the relationship of the structure $R$ to the energy $U$ of the system of interest. However, a potential energy function alone does not make a force field. It is the combination of the potential energy function with the parameters used in that function, as described below, that yield a force field;
%
\begin{eqnarray}
\label{potf}
U(R)=\sum_{bonds} K_b(b-b_0)^2 
+ \sum_{angles} K_{\theta}(\theta - \theta_0)^2 \nonumber\\
+ \sum_{dihedral} K_{\chi} \left(1 + \cos(n\chi - \delta)\right)  
+ \sum_{impropers} K_{imp}(\vartheta - \vartheta_0)^2 \nonumber\\
+ \sum_{nonbonded} \left(\epsilon_{ij}\left[\left(\frac{Rmin_{ij}}{r_{ij}}\right)^{12} - 2\left(\frac{Rmin_{ij}}{r_{ij}}\right)^6\right]\right) 
+ \frac{q_iq_j}{\epsilon r_{ij}}.
\end{eqnarray}
%
The equation \ref{potf} is an example of a potential energy function. It is comprised of a collection of simple functions to represent a minimal set of forces that can describe molecular structures. Bonds, angles, and out-of-plane distortions (improper dihedral angles) are treated harmonically, while dihedral or torsional rotations are described by a sinusoidal term. Interactions between atoms use a Lennard-Jones (LJ) 6-12 term to describe the atom-atom repulsion and dispersion interactions combined with electrostatics treated via a Coulombic term.\cite{Mackerell2004} 

In Eq.\ref{potf}, $b$ is the bond length; $\theta$ is the valence angle; $\chi$ is the dihedral or torsional angle; $\vartheta$ is the improper angle; $r_{ij}$ is the distance between atom $i$ and $j$. Parameters are the terms that represent the actual force field. Included in them are the bond force constant and the equilibrium distance, $K_b$ and $b_0$, respectively; the valence angle force constant and equilibrium angle, $K_{\theta}$, and $\theta_0$, respectively; the dihedral force constant, multiplicity and phase angle, $K_\chi$, $n$, and $\delta$, respectively; the improper force constant and equilibrium improper angle, $K_{imp}$ and $\vartheta_0$, respectively. Collectively, these represent the internal or intra-molecular parameters. Non-bonded parameters between atoms $i$ and $j$ include the partial atomic charges, $q_i$, $q_j$ and the LJ well-depth, $\epsilon_{ij}$, and minimum interaction radius, $Rmin_{ij}$, used to treat the van der Waals interactions. These terms are also referred to as interaction or external parameters. Typically, $\epsilon$ and $Rmin$ are obtained for individual atom types either by fitting to an experimental data or from quantum chemical calculations, which are then combined to yield $\epsilon_{ij}$ and $Rmin_{ij}$ for the interacting atoms via combining rules.\cite{D.2005}
\section{Periodic boundary conditions}
The correct treatment of boundaries and boundary effects is crucial to simulation methods. In order to simulate bulk phases it is essential to choose boundary conditions that mimic the presence of an infinite bulk surrounding the $N$ particle model system. This is achieved by employing periodic boundary conditions. The volume containing the $N$ particles is treated as the primitive cell of an infinite periodic lattice of identical cells. If a particle leaves the box during the simulation, then it is replaced by an image particle that enters from the opposite side. The number of particles within the central box thus remains constant. The intermolecular interactions which are of short range are truncated beyond a certain cutoff distance.
\section{Long-range interactions}
Intermolecular Lennard-Jones interactions are of short range and decay ($1/r^6$) rapidly with the distance $r$. Truncating these interactions to zero beyond a certain cutoff does not induce large errors. However, this is not the case for Coulombic interactions, which vary strongly at small distances and decay slowly at long distances. The Ewald summation is one of the methods by which long range interactions can be handled appropriately. The Ewald summation method splits the slowly converging sum over the Coulomb potential into two exponentially converging sums
\begin{equation}
\label{eq:ew1}
	\frac{1}{r}=\frac{f(r)}{r}+ \frac{1-f(r)}{r}.
 \end{equation}
The splitting function $f$ is chosen in a way that the first part $\frac{f(r)}{r}$ should be negligible or zero, beyond a certain cutoff distance. The second part $\frac{1-f(r)}{r}$ should be a slowly varying function of $r$, such that its Fourier transform can be represented by only a few $\mathbf{k}$ vectors. This permits an efficient calculation of this contribution in the reciprocal space. The function $f$ is a complementary error function 
\begin{equation}
\label{eq:ew2}
	erfc(r)=\frac{2}{\sqrt{\pi}}\int\limits^{\infty} _r \exp{(-t)^2} dt .
 \end{equation}
The potential energy of the primary box is now given by:
 \begin{eqnarray}
\label{eq:ew3}
V &=& V^r + V^k + V^s + V^d\\
V^r &=& \frac{1}{2}\sum\limits^{N}_{i=1}\sum\limits^{N}_{j=1}\sum_{n=0}{}^{'}\frac{q_{i}q_{j}}{4\pi\epsilon_0} \frac{erfc(\alpha\vert\mathbf r_{ij}+ \mathbf n\vert)} {\vert\mathbf r_{ij}+ \mathbf n\vert} \\
%
V^k &=& \frac{1}{2}\sum\limits_{k\neq 0}\sum\limits^{N}_{i=1}\sum^{N}_{j=1}\frac{1}{\pi L^3} \frac{q_{i}q_{j}}{4\pi\epsilon_0} \frac{4\pi^2}{k^2}\exp{\left( -\frac{k^2}{4\alpha^2}\right) cos(\mathbf{k} \cdot \mathbf{r}_{ij}) }\\
%
V^s &=& -\frac{\alpha}{\sqrt{\pi}}\sum^N_{k=1}\frac{q^2_k}{4\pi\epsilon_0} \\
%
V^d &=&\frac{2\pi}{3L^3}\left\vert\sum_{i=1}^N\frac{q_i}{4\pi\epsilon_0}\mathbf{r_i}\right\vert^2,
\end{eqnarray}
where $V^r$ is the contribution from real space, $V^k$ the contribution from the reciprocal space, $V^s$ the self energy and $V^d$ the dipole correction. The term $\alpha$ is often referred to as Ewald parameter and it tunes the relative weight of the real space and the reciprocal space contribution. The real space summation involving the error function converges very rapidly, and beyond a certain cutoff distance its value can be considered negligible. The sum of Gaussian functions in real space includes the interaction of each Gaussian with itself, thus it needs to be subtracted. In the reciprocal space summation, the vectors $\mathbf{k}$ are reciprocal vectors and are given by $\mathbf{k}=2\pi\mathbf{n}/L$. If the surrounding medium has an infinite relative permittivity, then no correction term is required. However, if the surrounding medium is vacuum, then $V^d$ must be added.\cite{Leach2001}
%
\section{Finite size effects}
Molecular simulations deal with a much smaller number of degrees of freedom ($ N \le 10^6$) than a typical experiment ($N \approx 10^{23}$), an aspect that can lead to artifacts in simulations. Finite system size effects are particularly large near phase transitions and percolation thresholds. In an infinite system, a first-order transition is characterized by a delta function singularity in the first derivatives of the thermodynamic potential, whereas for finite systems, delta function singularities are rounded off.
\par
The work of formation $\Delta G$ of a cluster of size $n$, according to classical nucleation theory, is given by:
%
\begin{equation}
\Delta{G} = -n\Delta\mu -\gamma{A},
\end{equation}
%
where $\Delta\mu$ is the difference in the chemical potentials of the vapor and the liquid phase, $\gamma$ is the surface tension and $A$ is the interfacial area. In a canonical ensemble with a finite number of particles far below the thermodynamic limit, the formation of a new phase causes a rapid consumption of the metastable phase. It can cause discrepancies between the results obtained from simulations and those of experiments. MacDowell \textit{et al.}\cite{Macdowell:04a} have shown that for a system having a volume $V$ and the total density $\rho$, the free energy cost $\Delta{A}_{hom}$ required to increase the density above the coexistence vapor density $\rho^c_g$, and the free energy cost $\Delta{A}_{inh}$ if the excess vapor condenses to form a droplet, can be given by:
%
\begin{eqnarray}
\Delta{A}_{hom} = \mu_c(\rho - \rho^c_g)V + \frac{1}{2}\frac{V}{\kappa^c_g}\left( \frac{\rho - \rho^c_g}{\rho^c_g}\right)^2 \\
\Delta{A}_{inh} = \mu_c(\rho - \rho^c_g)V + k_g\gamma\left(\frac{\rho - \rho^c_g}{\rho^c_l - \rho^c_g}\right)^{2/3}V^{2/3},
\end{eqnarray}
%
where $\mu_c$ and $\kappa^c_g$ are the chemical potential and the vapor compressibility at coexistence, $k_g$ is the geometric constant that allows the surface to be expressed in terms of $V^{2/3}$, and $\rho^l_c$ is  the coexistence liquid density. For small deviations from $\rho^c_g$, the equilibrium state of the system is a supersaturated homogeneous vapor as $\Delta{A}_{hom} < \Delta{A}_{inh}$. It must be emphasized that the system is in equilibrium and not in a metastable state. For larger deviations from $\rho^c_g$ , $\Delta{A}_{hom}$ becomes greater than $\Delta{A}_{inh}$ and condensation takes place. This is also a system size effect and it vanishes in the thermodynamic limit. It has also been shown that the density at which condensation takes place, $\rho^t$, scales as $\rho^t - \rho^c_g \propto V^{-1/4}$; hence, the transition moves towards coexistence with increasing the system size.
\par
%Most of the simulations of peptide aggregation are carried out in either a canonical or an  isothermic-isobaric ensemble. The number and concentration of the peptides should have a great impact on the stability of the peptide aggregate formed, which might be unstable due to less number of peptides available, even though the simulations were carried out at concentrations deep inside the co-existing region. However, the simulations of polymers in grand canonical ensembles are extremely difficult, as compared to the corresponding simulations of simple fluids. Hence, systematic studies on the effect of concentration and volume on peptide aggregation are extremely important.
\section{Methods}
%
\subsection{Molecular Dynamics simulations}
%
Molecular Dynamics simulations solve Newton's equation of motion of molecules to generate new configurations. For a system of $N$ interacting particles:
%
\begin{equation}
m_{i}\frac{\partial^2 \mathbf r_i}{\partial t^2} = \mathbf F_i,\ i=1 \cdots N,
\end{equation}
where $m_{i}$ is the mass of the particle $i$, $\mathbf F_i$ is the force and $\mathbf r_i$ is the position vector of the particle $i$. The forces are the negative derivative of given potential functions $U(\mathbf r_1, \mathbf r_2,\cdots,\mathbf r_N):$
\begin{equation}
\mathbf F_i= -\frac{\partial U}{\partial \mathbf r_i}.
\end{equation}
%
Solving Newton's equation of motion requires a numerical procedure for integrating the differential equation. Finite difference techniques are mainly used for this purpose. These techniques are used to generate molecular dynamics trajectories with continuous potential models, which is assumed to be pairwise additive. The integration is broken down into many small stages, each separated in time by a fixed time $\delta t$. The total force on each particle in the configuration at a time $t$ is calculated as the vector sum of its interactions with other particles. From the force, accelerations of particles are determined, which are then combined with the positions and velocities at a time $t$ to calculate the positions and velocities at a time $t + \delta t$. The force is assumed to be constant during the time step, which is generally in the range of 0.5--2.0 femtoseconds. The forces on the particles in their new positions are then determined, leading to new positions and velocities at time $t + 2\delta t$, and so on. The \emph{leap frog} modification of the \emph{Verlet} algorithm is probably the most widely used method for integrating the equations of motion. The velocities $\mathbf v(t+\frac{1}{2}\delta t)$ are first calculated from the velocities at time $t-\frac{1}{2}\delta t$ and acceleration at time $t$, as follows:
%
\begin{equation}
\mathbf v(t+\frac{1}{2}\delta t) = \mathbf v(t-\frac{1}{2}\delta t) +\delta t\mathbf a(t),
\end{equation}
the positions $\mathbf r(t + \delta t)$ are then deduced from velocities just calculated together with the positions at time $\mathbf r(t)$ using Equation \ref{leap}. Thus, the velocities \emph{leapfrog} over the position (hence the name)\cite{Allen1989} 
%
\begin{equation}
\label{leap}
\mathbf r(t +\delta t) = \mathbf r(t) + \delta t\mathbf v(t + \frac{1}{2}\delta t).
\end{equation}
%****************************************
\subsection{Monte Carlo simulations}
The configurational part of the partition function can also be denoted by $Z$:
\begin{equation}
\label{eq:mc1}
	Z \equiv \int d\mathbf r^N \exp[-\beta U(\mathbf r^N)],
\end{equation}
where $\beta=1/k_BT$ and $\mathbf r^N$ stands for coordinates of all $N$ particles. If random points are generated in the configurational space according to the probability density $\varrho(\mathbf r^N)$, where
\begin{equation}
\label{eq:mc2}
	 \varrho(\mathbf r^N)\approx \frac{\exp[-\beta U(\mathbf r^N)]}{Z}\ ,
\end{equation}
then the average property of a variable $A$ is:
\begin{equation}
\label{eq:mc3}
	 \langle A \rangle \approx \frac{1}{L}\sum\limits^{L}_{i=1}n_i A(\mathbf r^N),
\end{equation}
where the number of points $n_i$ generated per unit volume around a point $\mathbf r^N$ is equal to $L\varrho(\mathbf r^N)$ and $L$ is the total number of points generated.
The main idea of importance sampling is to generate points in configurational space with a relative probability proportional to the Boltzmann factor. If $o$ and $n$ are two configurations of the system, then the transition probability to go from $o$ to $n$, $\pi(o\rightarrow n)$ is:
\begin{equation}
\label{eq:mc4}
	 \pi(o\rightarrow n) = \alpha(o\rightarrow n) \times acc(o\rightarrow n),
 \end{equation}
where $\alpha(o\rightarrow n)$ is the probability to perform a trial move from $o$ to $n$ and $acc(o\rightarrow n)$ is the probability of accepting such a trial move. In equilibrium, the average number of accepted trial moves that results in the system leaving the state $o$ must be equal to the number of accepted trial moves from all other states $n$ to the state $o$. Generally, the detailed balance condition is applied:
\begin{equation}
\label{eq:mc5}
	 \varrho(o)\pi(o\rightarrow n) = \varrho(n)\pi(n\rightarrow o),
 \end{equation}
 where the average number of accepted moves from $o$ to $n$ is exactly canceled by the number of reverse moves.
 In the Metropolis scheme, the transition probability of going from $o$ to $n$ is given by:
 
 \parbox{8cm}{\begin{eqnarray*}
\label{eq:mc6}
	\pi(o\rightarrow n) & = & \alpha(o\rightarrow n)\\
	& = &\alpha(o\rightarrow n)[\varrho(n)/{\varrho(o)}] \nonumber 
\end{eqnarray*}}
\hfill
\parbox{4cm}{\begin{eqnarray}
	\varrho(n) \geq \varrho(o)\\ \nonumber
	\varrho(n) < \varrho(o).
\end{eqnarray}} \\
If $\alpha$ is chosen to be symmetric, then the equation \ref{eq:mc5} can be rewritten as:
 \begin{eqnarray}
 \label{eq:mc7}
 \varrho(o) acc(o\rightarrow n) &=& \varrho(n) acc(n\rightarrow o)\\ \nonumber
 \frac{acc(o\rightarrow n)}{acc(n\rightarrow o)} &=& \frac{\varrho(n)}{\varrho(o)} \\ \nonumber
  &=& \exp\{-\beta[U(n)-U(o)]\}.
 \end{eqnarray}
The metropolis scheme is implemented by attempting a trial displacement and calculating the change in energy $\Delta U=[U(n)-U(o)]$. If $\Delta U <0$ then the move is accepted. Otherwise a random number $X$ is generated and the move is accepted if $\exp(-\beta\Delta U) \ge X$.
%
\subsection{Replica Exchange simulations}
In order to carry out meaningful simulation studies of biomolecular systems containing thousands of atoms, a sampling of a huge configurational space is required, and therefore such studies are extremely time consuming. This problem is more acute at ambient and low temperatures, as the system could get trapped in local energy minima and might not be able to escape from it within the duration of the simulation run. Similar problems exist in other fields of physics and a parallel tempering technique (or the replica exchange algorithm) was proposed.\cite{Swendsen:86a,Marinari:92a} This method was applied to the simulation of biomolecules in 1997,\cite{Hansmann:97a} and extended to molecular dynamics simulations.\cite{Sugita1999} 
\par
Replica Exchange Molecular Dynamics (REMD) is an enhanced sampling technique based on the parallel tempering Monte Carlo method where multiple copies (or replicas) of identical systems are simulated in parallel at different temperatures.\cite{Sugita1999} Periodically state-exchange move are attempted, where two neighboring replicas exchange their thermodynamic states (their temperatures). The acceptance rule for each state-exchange move between two neighboring states $i$ and $j$ is chosen to be
\begin{equation}
\label{eq:rep1}
	 P(acc) = min\left\lbrace 1,\exp(1/k_BT_i -1/k_BT_j)(U(r_i^N)-U(r_j^N))\right\rbrace .
 \end{equation}
The state-exchange acceptance probability has been shown to obey the detailed balance condition for an extended ensemble of canonical states.\cite{smit2001}
