\documentclass{seminar}
\usepackage{fancybox}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\newcommand{\vect}[1]{\mbox{\boldmath $#1$}}

\begin{document}

%################################################################################
%################################################################################

\begin{slide*}
\pagestyle{empty}
\begin{center}
{\bf \large Modeling and simulating complex ionic systems in the presence of dielectric inhomogeneities}\\
\vspace{0.6cm}
\Ovalbox{{\it \large Dezs\H{o} Boda}}\\
\vspace{0.2cm}
{University of Veszprém, Department of Physical Chemistry}\\
\vspace{0.5cm}
\ovalbox{{\it Dirk Gillespie$^{1}$, Wolfgang Nonner$^{2}$, Douglas Henderson$^{3}$, Bob Eisenberg$^{1}$}}\\
\vspace{0.2cm}
{$^{1}$Rush University Medical Center, Chicago}\\
{$^{2}$University of Miami School of Medicine, Miami}\\
{$^{3}$Brigham Young University, Provo\\
\vspace{0.5cm}
December 13, 2004, Linz, Austria
}

\end{center}
\end{slide*}

%################################################################################


%################################################################################

\pagestyle{enyem}

\begin{slide*}
\slideheading{Molecular models - two particle system}
\parbox{4cm}{\framebox{\scalebox{0.35}{\includegraphics*{twobody.eps}}}\\
\begin{center}
\ovalbox{$u(\mathbf{r}_{12},\vect{\omega}_{1},\vect{\omega}_{2})$}
\end{center}
\begin{small}
\textbf{Intermolecular pair potential} 
\begin{itemize} 
  \item depends on the mutual position $(\mathbf{r}_{ij})$
  \item and the mutual orientations ($\vect{\omega}_{i}, \vect{\omega}_{j} $) of the molecules 
  \item It can be calculated from quantum mechanics or we can use \textbf{model potentials}.
\end{itemize}      
\end{small}}
\parbox{7.9cm}{
\raggedright
It is usual to \textbf{divide} the model potential into
\begin{itemize} 
  \item repulsive core (hard sphere, exponential, or $1/r^{12}$ repulsion)
  \item attractive dispersion force ($1/r^{6}$ or exponential Yukawa attraction)
  \item elecrostatic forces (Coulomb, dipole, quadrupole, etc.)
\end{itemize}      
\textbf{Example}: SPC/E  water
\begin{itemize}
\item a Lennard-Jones potential \ovalbox{$4\varepsilon\left[ \left( \dfrac{\sigma}{r}\right) ^{12} - \left( \dfrac{\sigma}{r}\right) ^{6}\right] $} on the oxygen
\item and partial charges on the atoms
\item 6 parameters: 2 LJ parameters, 2 partial charges, distance and angle of hydrogens
\end{itemize} 
}
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{Molecular models - many particle systems}
\parbox{4cm}{
\raggedright
\framebox{\scalebox{0.3}{\includegraphics*{manybody.eps}}}\\
\begin{small}
\textbf{Intermolecular potentials} acting between particles
\begin{itemize} 
  \item defines the system on the\textit{ microscopic level} 
  \item determines the \textbf{energy} for a given configuration of the particles
\end{itemize}      
\end{small}}
\parbox{7.9cm}{
\begin{itemize} 
  \item The energy is not equal to the sum of the pair potentials:
 	\ovalbox{$U \neq \sum_{i<j} u_{ij} $}
  \item The difference is due to \textbf{many body forces} 
  \item the vicinity of a third particle changes the pair potential
  \item \textbf{Solution:} using \textit{effective} pair potentials 
  \item fit the parameters of the potential to experimental results
  \item For instance: fit vapour-liquid equilibrium curves 
  \item \textbf{How can we relate the microscopic pair potential to mesaured physical quantities}?
\end{itemize}      

\begin{center}
\Ovalbox{\textbf{STATISTICAL MECHANICS}}
\end{center}

}
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{The relation of the microscopic and the macroscopic level}
\parbox{7cm}{
\raggedright
\framebox{\scalebox{0.4}{\includegraphics*{allen.eps}}}\vspace{0.3cm}
\begin{small}
\textbf{Computer simulations} serve as a bridge between real systems and their models
\end{small}}
\parbox{4.9cm}{
\raggedright
\begin{itemize} 
  \item \textbf{What is simulation? Theory? Experiment?}
  \item \textbf{Answer:} a statistical mechanical tool that gives \textit{exact results}
  \item \textbf{apart from} \textit{statistical uncertainty} and \textit{system size} problems.
  \item A \textbf{theory} has to be verified by comparison to simulation.
  \item A \textbf{model} has to be verified by comparing the simulation results to macroscopic experimental data.
\end{itemize}      
}
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{Statistical mechanics}
\raggedright
\begin{itemize} 
  \item The \textbf{thermodynamic state} has to be specified. For example: \textbf{canonical ensemble} - number of particles ($N$), volume ($V$), and temperature ($T$) are fixed.  
  \item The macroscopic quantities of the system are known if the partition function is known: \ovalbox{$Q=\sum_{i}\exp(-U_{i}/kT)$}, where we sum over the possible states of the system and $U_{i}$ is the energy in state $i$. In the classical theory the sum becomes an integral over phase space (position and velocity).
  \item For example: Helmholtz free energy: $F=-kT\ln Q$, etc. 
  \item \textbf{Ensemble average} of physical quantity $B$ is \ovalbox{$\langle B \rangle =\dfrac{1}{Q}\sum_{i}B_{i}\exp(-U_{i}/kT)$}, where $B_{i}$ is the value of $B$ in state $i$ - calculated by \textbf{Monte Carlo} (MC) simulation
  \item \textbf{Time average} of the physical quantity $B$ is  \ovalbox{$\bar{B}=\dfrac{1}{\tau}\int B(t)dt$} - \\calculated by \textbf{molecular dynamics} (MD) simulation
  \item \textbf{Postulate}: the two averages are equal - \ovalbox{$\langle B \rangle =\bar{B}$}.
\end{itemize}      
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{Molecular dynamics simulation}
\raggedright
\begin{itemize}
\item We move the particles according to the equations of motions of Newton - a deterministic approach.
\item The force acting on each particles is known if the intermolecular potentials are specified.
\item Choosing a small $\Delta t$ time interval, we can calculate the velocities and positions of the particles at time $t+\Delta t$ if we know the velocities and positions at time $t$.
\item In the simulation, we calculate and store the various physical quantities and we calculate the \textit{time average} at the end.
\item Appropriate to simulate \textbf{non-equilibrium} systems and \textbf{dynamical properties} (e. g. diffusion).
\item \textbf{Problems}: native ensemble is $ NVE $ (microcanonical).
\item Application to other ensembles: sort of mixing of MD and MC.
\item Problems with thermostating. 
\end{itemize}      
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{Monte Carlo simulations}
\raggedright
\textbf{A simplified picture}
\begin{itemize}
\item We randomly sample the phase space - hence the name
\item \textbf{First idea}: generate random possible states uniformly and estimate 
\begin{center}
$Q=\sum_{i}\exp(-U_{i}/kT)$
\end{center}
\item The ensemble average of $ B $ would be $\langle B \rangle =\sum_{i}B_{i}\mathsf{p}_{i}$.
\item where $\mathsf{p}_{i}=\dfrac{1}{Q}\exp(-U_{i}/kT)$ is the probability of state $ i $.
\item \textbf{Problem}: the number of possible states is huge.
\item We sample states that are far from equilibrium: $ U_{i} $ is large; $ \exp(-U_{i}/kT) $ is small.
\item Their contribution to $ Q $ is small - they are \textit{not important}.
\item \textbf{Solution}: \textit{importance sampling} - we sample states that are important
\end{itemize}      
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{Monte Carlo simulations}
\raggedright
\textbf{Importance sampling}
\begin{itemize}
\item Let us generate states with $ \mathsf{p}_{i} $ probability (Boltzmann sampling).
\item We have to bias in the ensemble average with the Boltzmann factor:\\
\begin{center}
\ovalbox{$ \langle B \rangle =\dfrac{\sum_{i}B_{i}\exp(-U_{i}/kT)/\mathsf{p}_{i}}{\sum_{i}\exp(-U_{i}/kT)/\mathsf{p}_{i}} =\dfrac{\sum_{i}B_{i}}{M} $}
\end{center} 
where $ M $ is the number of generated states.
\item This is a simple mean average.
\item Realization of this sampling: \textbf{Metropolis sampling} 
\\ \begin{small}N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, E. Teller, \textit{Chem. Phys.} 1953, \textbf{21}, 1087. (simulation of hard disks - Los Alamos)\end{small} 
\item There are a lot of other possibilities.	
\end{itemize}      
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{Monte Carlo simulations}
\raggedright
\textbf{Metropolis sampling}
\begin{itemize}
\item Probability of transition $i \rightarrow i+1 $ is 
\begin{center}
\ovalbox{$ \mathsf{p}_{i\rightarrow i+1} =\mathsf{p}_{i+1}/\mathsf{p}_{i}=\exp[-(U_{i+1}-U_{i})/kT] =\exp(-\Delta U/kT) $}
\end{center}
\item State $ i+1 $ is accepted by relating to the previous state $ i $ - a Markov chain.
\item The new state is accepted with the probability \ovalbox{$ \min[1,\exp(-\Delta U/kT)]$}.
\item The MC step is always accepted if $ \Delta U<0 $.
\item When $ \Delta U>0 $, it is accepted with probability $ \exp(-\Delta U/kT)  $
\item This is a Boltzmann sampling.
\end{itemize}      
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{Monte Carlo simulations}
\raggedright
\textbf{Other ensembles}
\begin{itemize}
\item $ NpT $: \textbf{fixed pressure} - volume change
\item $ \mu VT $ \textbf{grand canonical}: fixed chemical potential - particle insertion/deletion
\item \textbf{Gibbs ensemble} - two simulation boxes with particle and volume exchange - simulation of phase equilibria 
\item \textbf{Biased moves}: to improve the efficiency of sampling and accelerate simulations
\item \textbf{Example}: particle exchange between a small volume subsystem (ion channel) and a large volume subsystem (bath).
\item Normally, such moves are rare and sampling of occupancy of the channel is inefficient.
\end{itemize}      
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{Periodic boundary condition}
\parbox{7cm}{
\raggedright
\framebox{\scalebox{0.33}{\includegraphics*{pbc.eps}}}}
\parbox{4.9cm}{
\raggedright
\begin{itemize} 
  \item \textbf{imitation} of macroscopic (infinite) sample
  \item the Universe is filled with the \textbf{periodic images} of the central simulation cell
  \item if a particle moved out of the central simulation cell, its periodic image is coming back in from the other side
  \item \textbf{surface effects} are eliminated
  \item \textbf{other problems}: enforced periodicity, system size dependence prevails
\end{itemize}      
}
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{Mixed boundary conditions}
\parbox{7cm}{
\raggedright
\framebox{\scalebox{0.31}{\includegraphics*{dl.eps}}}
double layer (DL) geometry: PBC in $ x $ and $ y $ dimensions, confined by walls in $ z $ dimension
}
\parbox{4.9cm}{
\raggedright
\begin{itemize} 
  \item PBC is applied in dimensions where the system is \textbf{homogeneous}
  \item in the dimension of inhomogeneity the system can be confined by walls
  \item \textbf{Inhomogeneous systems}: the fluid at interfaces is in equilibrium with a bulk system.
  \item The bulk is a part of the simulation cell (far from walls) or specified by a chemical potential (grand canonical scheme). 
\end{itemize}      
}
\end{slide*}

%################################################################################
 
\begin{slide*}
\slideheading{Problems with large scale MD simulations}
\parbox{6.2cm}{
\raggedright
\framebox{\scalebox{0.22}{\includegraphics*{kcsa1.eps}}
					   \scalebox{0.20}{\includegraphics*{kcsa4.eps}}}
             \framebox{\scalebox{0.21}{\includegraphics*{kcsa2.eps}}
					   \scalebox{0.21}{\includegraphics*{kcsa3.eps}}}\\
\begin{small}D. Doyle \textit{et al.}, {\it Science} 1998, {\bf 280}, 69.\end{small}
}
\parbox{5.7cm}{
\raggedright
\begin{itemize} 
  \item If we know detailed 3D structure of the system (e. g. the KcsA potassium channel on the left), we can perform all atom MD simulation.
  \item LJ potentials and partial charges on atoms and ions, SPC/E for water.
  \item \textbf{For example:}\\ 
	\begin{small}T. W. Allen, S. Kuyucak, S. H. Chung: Molecular dynamics study of the KcsA potassium channel, \textit{Biophysical Journal}, 1999, \textbf{77}, 2502.\end{small}
 \item \textbf{Problems.} At least, in my view.  
\end{itemize}      
}
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{Problems with large scale MD simulations}
\parbox{6.2cm}{
\raggedright
\framebox{\scalebox{0.22}{\includegraphics*{kcsa1.eps}}
					   \scalebox{0.20}{\includegraphics*{kcsa4.eps}}}
             \framebox{\scalebox{0.21}{\includegraphics*{kcsa2.eps}}
					   \scalebox{0.21}{\includegraphics*{kcsa3.eps}}}\\
\begin{small}D. Doyle \textit{et al.}, {\it Science} 1998, {\bf 280}, 69.\end{small}
}
\parbox{5.7cm}{
\raggedright
\begin{itemize} 
   \item \textbf{Uncertainties of models and molecular parameters} - models and parameters that are good in a given situation might not be good in other situation.
  \item \textbf{Example}: water - some models are good for thermodynamics, others good for dielectric, interfacial, etc. properties.
  \item results should be verified by comparison to experimental data (current-voltage).
  \item long enough simulations cannot be run to get a macroscopic sample  
\end{itemize}      
}
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{Continuum models}
\raggedright
\begin{description}
\item[Summarized]: molecular simulations of systems modelled in full atomic detail \textbf{might be problematic:} \hspace{0.1cm}
      \begin{itemize} 
        \item \textit{time consuming} - many particles (macroscopic sampling is not possible)
        \item \textit{unknown 3D structures} (proteins)
        \item \textit{uncertainties in estimation of potential parameters}
      \end{itemize}      
\item["Coarse-graining"] \hspace{0.3cm}
      \begin{itemize} 
        \item {\it important} particles are treated explicitely (e. g. ions)
        \item {\it the rest} is smoothed into a continuum (e. g. solvent, proteins, membrane, electrode)
        \item the various components with different polarizabilities are modelled as continuums with different dielectric coefficients
	    \item it can be described by an inhomogeneous dielectric coefficient: $ \varepsilon(\textbf{r}) $ 
      \end{itemize}      
\item[Basic task:] solving the Poisson equation
\ovalbox{$-\varepsilon_{0}\nabla\cdot\left[  \varepsilon\left(  \mathbf{r}\right) \nabla
\psi\left(  \mathbf{r}\right)  \right]  =\rho\left(  \mathbf{r}\right)$}
\end{description}
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{Macroscopic electrostatics}
\raggedright
\begin{description}
\item[The total charge is split into two parts:] \hspace{0.3cm}
      \begin{itemize} 
        \item source (free or explicitely treated) charges: \ovalbox{$g(\mathbf{r})=\rho(\mathbf{r})/\varepsilon_{0}$} and
	    \item polarization charges \ovalbox{$h(\mathbf{r})=\rho_{pol}(\mathbf{r})/\varepsilon_{0} $} (\textbf{to be determined})
        \item constitutive relation: $\rho_{pol}(\mathbf{r})=-\nabla\cdot\mathbf{P}(\mathbf{r})$, where
        \item polarization: $ \mathbf{P}(\mathbf{r})=-\varepsilon_{0}\chi(\mathbf{r})\nabla\psi(\mathbf{r})$
        \item dielectric susceptibility: $\chi(\mathbf{r})=\varepsilon(\textbf{r})-1$
      \end{itemize} 
\item[The corresponding Poisson equation:] 
	$ \nabla^{2}\psi(\mathbf{r})=-\dfrac{1}{\varepsilon_{0}}[\rho(\mathbf{r})+\rho_{pol}(\mathbf{r})] $ \hspace{0.3cm}   
\item[The potential is also split into two parts:] \hspace{0.3cm}
	\begin{itemize} 
	\item 
\Ovalbox{$ \psi(\mathbf{r}) =\psi_{e}(\mathbf{r})+\psi_{i}(\mathbf{r}) = \int_{\mathcal{D}}G(\mathbf{r}-\mathbf{r}^{\prime})g(\mathbf{r}^{\prime})d\mathbf{r}^{\prime}+
\int_{\mathcal{D}}G(\mathbf{r}-\mathbf{r}^{\prime})h(\mathbf{r}^{\prime})d\mathbf{r}^{\prime} $}
	\item the Green-function satisfies $ \nabla^{2}G(\mathbf{r}-\mathbf{r}^{\prime})=-\delta(\mathbf{r}-\mathbf{r} ^{\prime}) $ 
	\item $ G( \mathbf{r}-\mathbf{r}^{\prime}) =1/ (4\pi |\mathbf{r}-\mathbf{r} ^{\prime}| )$ for the boundary condition $ \psi(\mathbf{r}) \rightarrow 0 $ if $ \mathbf{r}\rightarrow \infty $
	\end{itemize}
\end{description}
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{Variational formalism}
\raggedright
\begin{description}
\item[Solution of the Poisson equation minimizes the functional] \hspace{0.3cm}
\begin{itemize}
	\item \doublebox{$ I\left[  \psi\right]  =\dfrac{1}{2}\int_{\mathcal{D}}\nabla\psi\cdot\nabla\psi
d\mathbf{r}-\int_{\mathcal{D}}\psi\left[  g+\dfrac{1}{2}\nabla\cdot\left(
\chi\nabla\psi\right)  \right]  d\mathbf{r} $}
	\item In minimum the value $-\varepsilon_{0}I[\psi]$ is the electrostatic energy
 	\item Substituting the Green function form of $ \psi(\mathbf{r})$, the functional can be expressed in terms of the charges: $ I[g ,h] $
	\item At a fixed $g(\mathbf{r})$, we have to find $ h(\mathbf{r}) $ that minimizes $ I[g ,h] $.
	\item This is equivalent to minimizing the $ h $-dependent part of the functional $ I_{2}\left[  h\right]  $ for a fixed $ g $.
\end{itemize}
	\item [So, the task is to solve the extremum condition]
\doublebox{$\dfrac{\displaystyle \delta I_{2}[h]}{\displaystyle \delta h(\mathbf{r})}=0$}
\end{description}
\begin{itemize}
\item \begin{small}\underline{R. Allen, J.-P. Hansen, and S. Melchionna, \textit{Phys. Chem. Chem. Phys.} \textbf{3}, 4177 (2001).}\end{small}
\item Allen \textit{et al.} solved it with the \textit{steepest descent} method.
\end{itemize}
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{Solution of the extremum condition}
\raggedright
\begin{center}
\doublebox{$\dfrac{\displaystyle \delta I_{2}[h]}{\displaystyle \delta h(\mathbf{r})}=0$}
\end{center}
We present a different solution resulting in a matrix equation.
\textbf{The following integral equation can be obtained from the extremum condition:}
\begin{center}
\shadowbox{
$
\begin{array}{rl}
&  h(\mathbf{r})\varepsilon(\mathbf{r})-\int_{\mathcal{D}}h(\mathbf{r}%
^{\prime})\nabla_{\mathbf{r}}\varepsilon(\mathbf{r})\cdot\nabla_{\mathbf{r}%
}G(\mathbf{r}-\mathbf{r}^{\prime})d\mathbf{r}^{\prime}\\
&  =\nabla_{\mathbf{r}}\varepsilon(\mathbf{r})\cdot\nabla_{\mathbf{r}}\psi
_{e}(\mathbf{r})-\dfrac{1}{\varepsilon_{0}}\left[  \varepsilon(\mathbf{r}%
)-1\right]  \rho(\mathbf{r})
\end{array}
$
}
\end{center}
\begin{itemize} 
\item This equation should be discretized and solved for $ h $.
\item It is \textbf{general} and valid for arbitrary $ \varepsilon(\mathbf{r}) $ and $ \rho(\mathbf{r}) $.
\item We can make our life easier by making some assumptions for $ \varepsilon(\mathbf{r}) $ and $ \rho(\mathbf{r}) $.
\end{itemize} 
%\ovalbox{$ \psi_{e}(\textbf{r})=? $}
\begin{small}\underline{D. Boda, D. Gillespie, W. Nonner, D. Henderson, B. Eisenberg, \textit{Phys. Rev. E} 2004, \textbf{69}, 046702.}\end{small}
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{Point charges as source charges}
\raggedright
\begin{itemize} 
\item Suppose that the source charges are point charges: \ovalbox{$ \ovalbox{$\rho(\mathbf{r})=\sum_{k}z_{k}e\delta(\mathbf{r}-\mathbf{r}_{k})	 $} $}
\item Polarization charges of magnitude \ovalbox{$ -z_{k}e\dfrac{\varepsilon(\mathbf{r}_{k})-1}{\varepsilon(\mathbf{r}_{k})} \delta(\mathbf{r}-\mathbf{r}_{k}) $} are induced around the point charges (localized on them).
\item The electric potential produced by the point charges \textbf{and} the polarization charges induced on them: 
\ovalbox{$\psi_{e}(\mathbf{r})=\dfrac{\displaystyle e}{\displaystyle 4\pi\varepsilon_{0}}\sum_{k}\dfrac{\displaystyle z_{k}} {\displaystyle \varepsilon(\mathbf{r}_{k})|\mathbf{r}-\mathbf{r}_{k}|}$}\\ 
where $ \varepsilon(\mathbf{r}_{k}) $ is the dielectric coefficient at the position of the $ k $th source charge.
\item From now on, $ h(\textbf{r}) $  denotes \textbf{only} the polarization charges induced on the dielectric interfaces. In this case, the basic equation reduces to :\\

\shadowbox{
$
h(\mathbf{r})\varepsilon(\mathbf{r})-\int_{\mathcal{D}}h(\mathbf{r}^{\prime
})\nabla_{\mathbf{r}}\varepsilon(\mathbf{r})\cdot\nabla_{\mathbf{r}%
}G(\mathbf{r}-\mathbf{r}^{\prime})d\mathbf{r}^{\prime}=\nabla_{\mathbf{r}%
}\varepsilon(\mathbf{r})\cdot\nabla_{\mathbf{r}}\psi_{e}(\mathbf{r})$
}
\item This is valid for arbitrary $ \varepsilon(\mathbf{r}) $.
\end{itemize}
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{Sharp dielectric interfaces}
\raggedright
%\begin{center}
%$ \Longrightarrow $
%\ovalbox{$
%h(\mathbf{r})\varepsilon(\mathbf{r})-\int_{\mathcal{D}}h(\mathbf{r}^{\prime
%})\nabla_{\mathbf{r}}\varepsilon(\mathbf{r})\cdot\nabla_{\mathbf{r}%
%}G(\mathbf{r}-\mathbf{r}^{\prime})d\mathbf{r}^{\prime}=\nabla_{\mathbf{r}%
%}\varepsilon(\mathbf{r})\cdot\nabla_{\mathbf{r}}\psi_{e}(\mathbf{r})$
%}
%\end{center}
\begin{itemize} 
\item When the dielectric is piecewise uniform, $ \varepsilon(\mathbf{r}) $ jumps from one value to another along a boundary.
\item In this case, the volume integrals over $ \mathcal{D} $ become surface integrals over the dielectric boundary surfaces $ \mathcal{B} $.
\item The basic equation then reduces to:
\begin{center}
\shadowbox{$h(\mathbf{s})\varepsilon(\mathbf{s})-\Delta\varepsilon(\mathbf{s}%
)\int_{\mathcal{B}}h(\mathbf{s}^{\prime})\nabla_{\mathbf{s}}G(\mathbf{s}%
-\mathbf{s}^{\prime})\cdot\mathbf{n}(\mathbf{s})d\mathbf{s}^{\prime}%
=\Delta\varepsilon(\mathbf{s})\nabla\psi_{e}(\mathbf{s})\cdot\mathbf{n}%
(\mathbf{s})$}
\end{center}
where
\item $ h(\mathbf{s}) $ is now surface charge distribution
\item $ \varepsilon(\mathbf{s}) $ is the mean average of the dielectric constants on the two sides of the surface
\item $ \Delta\varepsilon(\mathbf{s}) $ is their difference, and 
\item $ \mathbf{n}(\mathbf{s})	 $ is a normal vector of the surface


\end{itemize}
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{The matrix equation}
\raggedright
%\begin{center}
%$ \Longrightarrow $\ovalbox{
%$h(\mathbf{s})\varepsilon(\mathbf{s})-\Delta\varepsilon(\mathbf{s}%
%)\int_{\mathcal{B}}h(\mathbf{s}^{\prime})\nabla_{\mathbf{s}}G(\mathbf{s}%
%-\mathbf{s}^{\prime})\cdot\mathbf{n}(\mathbf{s})d\mathbf{s}^{\prime}%
%=\Delta\varepsilon(\mathbf{s})\nabla\psi_{e}(\mathbf{s})\cdot\mathbf{n}%
%(\mathbf{s})$}
%\end{center}
The previous equation is valid for any discrete value of $ \mathbf{s}_{\alpha} $. Use the notations
\begin{itemize} 
\item normal vector: $\mathbf{n}_{\alpha}=\mathbf{n}(\mathbf{s}_{\alpha})	$
\item mean dielectric constant: $\varepsilon_{\alpha}=\varepsilon(\mathbf{s}_{\alpha})	$
\item jump in the dielectric constant: $\Delta\varepsilon_{\alpha}=\Delta\varepsilon(\mathbf{s}_{\alpha})$
\end{itemize}
The surface $ \mathcal{B} $ is discretized into surface elements. Let us suppose that the induced surface charge is constant on surface element $ \beta $ denoted by $ h_{\beta} $. Then our equation becomes
\begin{center}
\shadowbox{$\sum_{\beta}  h_{\beta} \left[ \varepsilon_{\beta} \delta_{\alpha\beta} - \Delta\varepsilon_{\alpha}  \int_{a_{\beta}}  \nabla_{\textbf{s}_{\alpha}} G(\textbf{s}_{\alpha}-\textbf{s}_{\beta}) \cdot \textbf{n}_{\alpha} d\textbf{s}_{\beta} \right]
= \Delta\varepsilon_{\alpha} \nabla \psi_{e} (\textbf{s}_{\alpha}) \cdot \textbf{n}_{\alpha} .
\label{surface2}$}
\end{center}
where the integral expresses the polarization of the surface element $ \beta $ by the induced charge at $ \mathbf{s}_{\alpha} $, and vice versa. The integral is taken for the $ \beta $th surface element.

%The surface integral becomes a sum, and we obtain:\\
%\begin{center}
%\shadowbox{$\sum_{\beta}h_{\beta}\left(  \delta_{\alpha\beta}\varepsilon_{\beta}%
%-\Delta\varepsilon_{\alpha}\frac{\displaystyle \partial G_{\alpha\beta}}{\displaystyle \partial n_{\alpha
%}}a_{\beta}\right)  =\Delta\varepsilon_{\alpha}\frac{\displaystyle \partial\psi_{\alpha}%
%}{\displaystyle \partial n_{\alpha}}\label{h_eq}$}
%\end{center}

\end{slide*}

%\begin{itemize} 
%\item vector of center of mass: $ \mathbf{s}_{\beta} $
%\item area: $a_{\beta}$
%\end{itemize}


%################################################################################

\begin{slide*}
\slideheading{The matrix equation}
\raggedright
%\begin{center}
%$ \Longrightarrow $\ovalbox{$\sum_{\beta}h_{\beta}\left(  \delta_{\alpha\beta}\varepsilon_{\beta}%
%-\Delta\varepsilon_{\alpha}\frac{\displaystyle \partial G_{\alpha\beta}}{\displaystyle \partial n_{\alpha
%}}a_{\beta}\right)  =\Delta\varepsilon_{\alpha}\frac{\displaystyle \partial\psi_{\alpha}%
%}{\displaystyle \partial n_{\alpha}}\label{h_eq}$}
%\end{center}
If $\mathbf{s}_{\alpha}$ is the center of the $ \alpha $th surface element, the previous equation can be written in the matrix form: \ovalbox{${\Huge A\mathbf{h}=\mathbf{c}} $}, where
\begin{itemize}
\item the \ovalbox{$A_{\alpha\beta}=\varepsilon_{\beta} \delta_{\alpha\beta} - \Delta\varepsilon_{\alpha}  \int_{a_{\beta}}  \nabla_{\textbf{s}_{\alpha}} G(\textbf{s}_{\alpha}-\textbf{s}_{\beta}) \cdot \textbf{n}_{\alpha} d\textbf{s}_{\beta}$} matrix depends only on the geometry of the dielectric
\item $ \mathbf{h} $ is the vector of induced charges 
\item \ovalbox{$c_{\alpha}=\Delta\varepsilon_{\alpha} \nabla \psi_{e} (\textbf{s}_{\alpha}) \cdot \textbf{n}_{\alpha}$} depends on the source charges
\item If the geometry does not change in a simulation, the matrix need be inverted only once at the beginning of the simulation. 
\item Moving the source charges (the ions) changes the vector $ \mathbf{c} $, from which the induced charges are calculated by the matrix-vector multiplication $ \mathbf{h}=A^{-1}\mathbf{c}  $.
\end{itemize}

\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{Calculation of the integral}
\raggedright
Regarding the calculation of the integral $ I_{\alpha\beta}=\int_{a_{\beta}}  \nabla_{\textbf{s}_{\alpha}} G(\textbf{s}_{\alpha}-\textbf{s}_{\beta}) \cdot \textbf{n}_{\alpha} d\textbf{s}_{\beta}$ there are two approaches 
\begin{enumerate}
\item Treating the surface charge $ h_{\beta} $ as a \textbf{point charge }of magnitude $ h_{\beta}a_{\beta} $ at the center of the $ \beta $th surface element of area $ a_{\beta} $
\begin{center}
\ovalbox{$I_{\alpha\beta} = \nabla_{\textbf{s}_{\alpha}} G(\textbf{s}_{\alpha}-\textbf{s}_{\beta}) \cdot \textbf{n}_{\alpha} a_{\beta}$}
\end{center}
for $ \beta \neq \alpha $ and $ 0 $ otherwise [\underline{Boda \textit{et al.}, \textit{Phys. Rev. E}, 2004, \textbf{69}, 046702.}].
This matrix has been developed before from the boundary condition $ D_{1n}=D_{2n} $, $ E_{1t}=E_{2t} $. [\underline{Hoshi \textit{et al.} JCP, 1987; Lu and Green, PCPS, 1997}]
\item Treating  $ h_{\beta} $ as a constant \textbf{surface charge} over the $ \beta $th surface element, and calculating $ I_{\alpha\beta} $ numerically after an appropriate parametrization of the surface [\underline{Nonner and Gillespie, \textit{Biophys. J.} in preparation}]. This approach is especially important if the surface is curved, therefore, we use the name \textit{curvature correction}.

\end{enumerate}

\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{Calculation of the energy}
\raggedright
\begin{itemize} 
\item The \textbf{source charge - source charge} interaction energy:
\ovalbox{$ W_{e}=\dfrac{1}{2}\sum_{j}ez_{j}\psi_{e}(\mathbf{r}_{j})$}
\item The \textbf{source charge - induced charge} interaction energy:
\begin{center}
\ovalbox{$ W_{i} = \dfrac{e}{8\pi}  \sum_{j} z_{j} \sum_{\beta} h_{\beta} \int_{a_{\beta}} \dfrac{d\mathbf{s}}{|\mathbf{s}-\mathbf{r}_{j}|}$}
\end{center}
Similarly to the calculation of the integral $ I_{\alpha\beta} $, this energy can also be calculated on two levels of approximation:
\begin{enumerate}
\item Treating the induced charge as a\textbf{ point charge }of magnitude $ h_{\beta}a_{\beta} $, the integral in the above equation becomes $ \dfrac{a_{\beta}}{|\mathbf{s}_{\beta}-\mathbf{r}_{j}|} $.
\item Treating the induced charge as a \textbf{surface charge} with the constant value $ h_{\beta} $ over the $ \beta $th surface element, the integral should be calculated numerically. This is a \underline{time consuming step in a simulation} because the $  |\mathbf{s}-\mathbf{r}_{j}|$ distances have to be evaluated once the $ j $th ion is displaced.
\end{enumerate}

\end{itemize}

\end{slide*}


%################################################################################

\begin{slide*}
\slideheading{Ion is modelled as a hard sphere with a point charge in the center} 
\raggedright
\begin{itemize} 
\begin{small}\item The dielectric constant inside the hard spheres is the same as outside.
\item The ions cannot overlap with the dielectric boundaries and cannot leave their host dielectric.
\item Otherwise, new dielectric boundaries would appear - matrix $ A $ should be inverted again.\end{small}
\end{itemize}
\begin{center}
\scalebox{0.35}{\includegraphics*{ionmodel.eps}}
\end{center}
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{One flat interface: $\varepsilon_{1} | \varepsilon_{2}$ geometry - ions on the right side}
\raggedright
\begin{center}
Comparison of density profiles obtained from the ICC and the image charge method \\
for various resolutions of the grid. [\underline{Boda \textit{et al.} \textit{PRE}, 2004, \textbf{69}, 046702.}]\\
\vspace{0.5cm}
\scalebox{0.4}{\includegraphics*{flat1.eps}}
\scalebox{0.4}{\includegraphics*{flat2.eps}}
\end{center}
\end{slide*}
%################################################################################

\begin{slide*}
\slideheading{Two flat interfaces: $\varepsilon_{1} | \varepsilon_{2} | \varepsilon_{3}$ geometry - electrical double layer}
\raggedright
It is usual to assume that the inner layer has a dielectric constant different from that of the diffuse layer or the elecrode because of the need for empirical adjustments to obtain agreement with experiments.
%By choosing an appropriate value for $ \varepsilon_{2} $, the capacitance of the double layer can be fitted to experimental values.
The elecrostatic consistency of the approach should be studied.
\begin{center}
%\vspace{0.5cm}
\scalebox{0.38}{\includegraphics*{dlgeometry.eps}} \hspace{0.5cm}
\scalebox{0.34}{\includegraphics*{dlprofile.eps}}
\end{center}
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{A dielectric sphere $\varepsilon_{1}=80$ embedded in a dielectric $\varepsilon_{2}=2$}
\raggedright
\begin{itemize} 
\item \underline{Radius of the sphere}: $R=5$ {\AA}. 
\item \underline{Distance of a source charge $ e $ from the center of the sphere}: $d=4$ {\AA}
\item \underline{Grid}: spherical coordinates $\vartheta$ and $\phi$ evenly divided into intervals $\Delta\theta = \Delta\phi$, $ N $ is the number of surface elements.
\item \underline{Analytic solution}: series in terms of Legendre polynomials.
\end{itemize}

%Ez 32 indukált töltést jelent.\\

\begin{center}
\scalebox{0.3}{\includegraphics*{1sphere.eps}} \hspace{1cm}
\scalebox{0.3}{\includegraphics*{sphere-grid.eps}}
\end{center}	
\end{slide*}

%################################################################################

\begin{slide*}
%\slideheading{Charge inside a dielectric sphere}
%Dielectric sphere $\varepsilon_{1}=80$ embedded in a dielectric $\varepsilon_{2}=2$
\begin{center}
\scalebox{0.43}{\includegraphics*{provo1.eps}}
\end{center}
The \textit{curvature correction} is important!!! 
\end{slide*}
%################################################################################

\begin{slide*}
\slideheading{The effective force acting between two charged dielectric spheres}
\raggedright
The effective pair potential is defined as \ovalbox{$ \varphi(r) =W(r)-W(r\rightarrow\infty) $}. The energy $ W(r\rightarrow\infty) $ corresponds to the interaction energy between the charge and the induced charge on its own sphere. The effective force is \ovalbox{$ F(r)=-d\varphi(r)/dr $}. 
\begin{center}
\scalebox{0.25}{\includegraphics*{fig6.eps}} \hspace{0.3cm}
\scalebox{0.43}{\includegraphics*{fig7.eps}}
\end{center}
%The interaction between the spheres becomes weaker as the dielectric constant inside the spheres decreases. 
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{MC simulation for a model calcium channel}
\raggedright
The selectivity filter of the calcium channel is known to be lined by 4 glutamate groups (EEEE locus). They are represented by 8 half charged oxygen ions that are confined to the filter but can move freely inside. A 0.1M NaCl solution used and a few Ca$ ^{++} $ ions. \textbf{Question:} whether the Na$ ^{+} $ or the Ca$ ^{++} $ enters the filter with a higher probability?
\begin{center}
\rotatebox{0}{\scalebox{0.40}{\includegraphics*{geometry.eps}}} \hspace{0.2cm}
\frame{\scalebox{0.2}{\includegraphics*{grid.eps}}}
\end{center}
\end{slide*}

%################################################################################


\begin{slide*}
\slideheading{MC simulation for a model calcium channel}
\parbox{6.2cm}{
\raggedright
\framebox{\scalebox{0.37}{\includegraphics*{calcium-selectivity.eps}}}
Results for uniform dielectric constant - $ \varepsilon =78.5 $ everywhere, no vestibule. \\
\begin{small}\underline{D. Boda \textit{et al.}, {\it Molec. Phys.}, 2002 {\bf 100}, 2361.}\end{small}
}
\parbox{5.7cm}{
\raggedright
\begin{itemize} 
\item The concentration of sodium ions is kept constant in the bath at 0.1 M. CaCl$ _{2} $ is added to the system.
\item The calcium ions replace sodium ions even when their bath concentration is 0.01 times smaller than that of the sodium ions.  
\item \textbf{Explanation}: \textit{charge-space competition} (CSC) mechanism  
\item The calcium ions provide two times the charge to neutralize the glutamate groups than sodium ions using about the same volume.
\end{itemize}      
}
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{Preliminary results for inhomogeneous dielectric}
%We have about 100 times more Na$ ^{+} $ ion in the bath than Ca$ ^{++} $ ion.
\begin{center}
\scalebox{0.41}{\includegraphics*{lg1.9.eps}}
\end{center}
\end{slide*}

%################################################################################

\begin{slide*}
\begin{center}
\frame{\scalebox{0.22}{\includegraphics*{CaNa.eps}}}
\end{center}
\end{slide*}

%################################################################################

\begin{slide*}
\begin{center}
\frame{\scalebox{0.22}{\includegraphics*{1Ca-b.eps}}}
\end{center}
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{Problem of absence of charge neutral bulk}
\parbox{6.2cm}{
\raggedright
\framebox{\scalebox{0.34}{\includegraphics*{no-bulk.eps}}}
The previous figure with a different scale showing the bulk. There is \textbf{no} an electroneutral bulk.
}
\parbox{5.7cm}{
\raggedright
\textbf{Explanation:}
\begin{itemize} 
\item \textbf{``Fixed charge ensemble''}
\item Normal procedure in \textbf{homogeneous} $ \varepsilon $ simulations: the number of ions is chosen so that the system is electroneutral.
\item The electroneutral bulk is automatically formed far from the walls (inhomogeneities).
\item Here, polarization charges are induced on the dielectric boundaries and \textbf{prevent} the formation of charge neutral bulk.
\end{itemize}      
}
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{Simulation by imposing Dirichlet or Neumann boundary conditions}
\parbox{6.7cm}{
\raggedright
\framebox{\scalebox{0.34}{\includegraphics*{newcell.eps}}}
Idea of Nonner and Gillespie (still in preparation): 
Define confining surface of the system (they call it electrode).
Apply an appropriate surface charge $ \sigma (\mathbf{r}) $ on the electrode to fulfill BC.
}
\parbox{5.2cm}{
\raggedright
\begin{itemize} 
\item \textbf{Dirichlet} BC: prescribed potential on the surface: $ V(\mathbf{s}) $
\item \textbf{Neumann} BC: prescribed normal electric field on the surface: $ E_{n} (\mathbf{s}) $
\item As the ions move, both the polarization charge $ h(\mathbf{s}) $ and the electrode charge $ \sigma(\mathbf{s}) $ are changing.
\item Electrode charge induces polarization charge and vice versa, so their calculation has to be coupled.
\end{itemize}      
}
\end{slide*}

%################################################################################

\begin{slide*}
\slideheading{Coupling induced and electrode charge calculations}

Equation for the dielectric boundaries $ \mathcal{B} $:

\begin{center}
\Ovalbox{$
\begin{array}{rcl}
4\pi\dfrac{\varepsilon(\textbf{s})}{\Delta\varepsilon(\textbf{s})}h(\textbf{s}) 
&+& \int_{ \mathcal{B}} h(\textbf{s}^{\prime}) \dfrac{(\textbf{s}-\textbf{s}^{\prime})\cdot \textbf{n}(\textbf{s})}{|\textbf{s}-\textbf{s}^{\prime}|^{3}}d\textbf{s}^{\prime}
+ \int_{ \mathcal{S}} \sigma(\textbf{s}^{\prime}) \dfrac{(\textbf{s}-\textbf{s}^{\prime})\cdot \textbf{n}(\textbf{s})}{|\textbf{s}-\textbf{s}^{\prime}|^{3}}d\textbf{s}^{\prime} \\
&=& -\sum_{k} q_{k}\dfrac{(\textbf{s}-\textbf{r}_{k})\cdot \textbf{n}(\textbf{s}) }{\varepsilon(\textbf{r}_{k})|\textbf{s}-\textbf{r}_{k}|^{3}}
\end{array}
$}
\end{center}
\vspace{0.5cm}
Equation for the electrode $ \mathcal{S} $:

\begin{center}
\Ovalbox{$\int_{ \mathcal{B}} \dfrac{h(\textbf{s}^{\prime})}{|\textbf{s}-\textbf{s}^{\prime}|} d\textbf{s}^{\prime} 
+\ \int_{\mathcal{S}} \dfrac{\sigma(\textbf{s}^{\prime})}{|\textbf{s}-\textbf{s}^{\prime}|}  d\textbf{s}^{\prime} 
=4\pi\varepsilon_{0}V(\textbf{s})- 
\sum_{k} \dfrac{q_{k}}{\varepsilon(\textbf{r}_{k})|\textbf{s}-\textbf{r}_{k}|}$}
\end{center}
\vspace{0.5cm}
To solve the problem, both surfaces $ \mathcal{B} $ and $ \mathcal{S} $ have to be discretized resulting in a matrix equation. \textit{Work is in progress...} 
\end{slide*}
%################################################################################

\begin{slide*}
\slideheading{Acknowledgement}
\begin{itemize} 
\item Dr. M\'onika Valisk\'o (University of Veszpr\'em)
\item Dr. David Busath (Brigham Young University)
\item T\'{\i}mea Nagy (University of Veszpr\'em)
\item Tibor Varga  (University of Veszpr\'em)
\end{itemize}
\vspace{0.5cm}
\begin{center}
\textbf{Download:} http://egsys.hu/$ ^{\sim} $dezo/Citations/citations.html
\end{center}
\vspace{0.5cm}
\begin{center}
\begin{large}\textbf{Thank you for your attention!}\end{large}
\end{center}
\end{slide*}

\end{document}


