1521 lines
137 KiB
TeX
1521 lines
137 KiB
TeX
|
||
% this file is called up by thesis.tex
|
||
% content in this file will be fed into the main document
|
||
|
||
%: ----------------------- introduction file header
|
||
%-----------------------
|
||
%%\usepackage[tbtags]{amsmath}
|
||
\newcommand{\boldm}[1] {\mathversion{bold}#1\mathversion{normal}}
|
||
%\usepackage{mathrsfs}
|
||
%
|
||
%\usepackage{bm}
|
||
% the code below specifies where the figures are stored
|
||
\ifpdf
|
||
\graphicspath{{2_Introduction/figures/PNG/}{2_Introduction/figures/PDF/}{2_Introduction/figures/}}
|
||
\else
|
||
\graphicspath{{2_Introduction/figures/EPS/}{2_Introduction/figures/}}
|
||
\fi
|
||
|
||
%----------------------------------------------------------------------
|
||
%: ----------------------- introduction content -----------------------
|
||
%----------------------------------------------------------------------
|
||
|
||
%: ----------------------- HELP: latex document organisation
|
||
% the commands below help you to subdivide and organise your thesis
|
||
% \chapter{} = level 1, top level
|
||
% \section{} = level 2
|
||
% \subsection{} = level 3
|
||
% \subsubsection{} = level 4
|
||
% note that everything after the percentage sign is hidden from output
|
||
|
||
\chapter{Computational Methods} \label{chap:comput_method}
|
||
|
||
Chemistry is a central scientific discipline dealing with the construction, transformation and properties of molecules. It provides a foundation for
|
||
understanding both basic and applied scientific disciplines at a fundamental level.\cite{Brown2009} Chemistry started out as an experimental
|
||
science. So, historically, all the chemical processes studied and the theories developed have been through experience. However, with the
|
||
development of theoretical concepts dealing with the structure of matter and with recent developments in computer technology, a new trend
|
||
in chemical research has developed: \textbf{theoretical chemistry}. \textbf{Computational chemistry}, as a subfield of theoretical chemistry,
|
||
combines fundamental laws with mathematical methods to solve chemical problems through calculations or simulations. Nowadays, modelling
|
||
of physico-chemical properties has become a important part of research in science of matter and the role of computational chemistry has become
|
||
fundamental. However, when resorting to computational chemistry, one always has to choose \textbf{a balance between accuracy and computational
|
||
efficiency}.
|
||
|
||
There exist many methods to model the interactions between different particles in atomic or molecular aggregates. Some of them are depicted in
|
||
Figure~\ref{methods} and are classified in terms of system sizes (\textit{y}-axis) and simulation times (\textit{x}-axis) they can tackle. If focusing
|
||
on the \textit{y}-axis, the closer a given method is from the origin of the graph, the more accurate it is. The farther from the origin a method is,
|
||
the larger are the system sizes it can simulate. For instance, \textbf{force fields} (FF), also referred to as \textbf{molecular mechanics} (MM)
|
||
methods, describe the interactions between particles by empirical interatomic potentials and electrons are not treated explicitly. FF methods thus
|
||
bypass the solving of the \textbf{Schr{\"o}dinger equation} and the quantum aspects of nuclear motion are neglected which allow them to model
|
||
very large systems. On the opposite, electronic structures methods such as \textbf{full configuration interaction} (Full CI) make it possible to
|
||
describe electron distributions with a very high accuracy but only for very small systems.
|
||
|
||
\begin{figure}[h]
|
||
\begin{center}
|
||
\includegraphics[width=12cm]{methods.png}
|
||
\caption{Comparison of the computational efficiency, \textit{i.e.} system sizes and simulation times, of various computational chemistry methods.
|
||
The \textit{y}-axis indicates the length of time accessible from classical molecular simulations for average system sizes tackle by each
|
||
method. The \textit{x}-axis indicates the approximative maximum system size tractable by each method in a single-point energy
|
||
calculation.} \label{methods}
|
||
\end{center}
|
||
\end{figure}
|
||
|
||
In parallel to the accuracy, the ability of a given method to efficiently and accurately explore PES is also of paramount importance
|
||
to model a variety of properties. Here also, different methods allow for different possibilities. For instance, the multi-configuration time-dependent Hartree
|
||
(MCTDH) method allows for a full quantum treatment of the nuclear degrees of freedom as long as only a limited number of them is taken into account.\cite{MCTDH}
|
||
For large systems, to access time-dependent information, one has to resort to \textbf{molecular dynamics simulations} which is dealt with \textbf{classical mechanics},
|
||
\textit{i.e.} Newton's second law.
|
||
|
||
This chapter gives a brief description of the theoretical foundations of the computational methods that were used during this thesis. Two main aspects are
|
||
developed. First, I describe the main methods to compute PES, \textit{i.e.} the electronic energy, through the solving of the Schr{\"o}dinger equation. Second,
|
||
I present the main computational tools that I have used to explore PES.
|
||
|
||
\section{Schr{\"o}dinger Equation}
|
||
|
||
The discovery of the Schr{\"o}dinger equation by E. Schr{\"o}dinger in 1925 is an extremely significant landmark in the development of the quantum mechanics. The behavior of a molecular system can be described by the Schr{\"o}dinger equation (a linear partial differential equation), which describes the wavefunction or state function of a quantum-mechanical system.\cite{Griffiths2018}
|
||
In quantum mechanics, the concept of wavefunction $\Psi$ is a fundamental postulate which defines the state of a system at each spatial position and time. %With these postulates, Schr{\"o}dinger equation can be derived from the fact that the time-evolution operator must be unitary and be generated by the exponential of a self-adjoint operator which is equivalent to a Hermitian operator $\hat{H}$ in the finite-dimensional case.
|
||
The Schr{\"o}dinger equation governs the evolution of the wavefunction $\Psi$ of particles in an atomic or molecular system.
|
||
%The form of the Schr{\"o}dinger equation depends on the physical situations.
|
||
In the case of a system composed of $M$ nuclei and $N$ electrons, the Schr{\"o}dinger equation for the wavefunction in position space $\Psi(\textbf{R$_\alpha$}, \textbf{r$_j$}, t)$ can be written as:
|
||
%\begin{eqnarray}
|
||
%\label{TDSE}
|
||
%&i\hbar\frac{\partial\Psi(\textbf{\boldm{$R_1$}}, \mathbf{R_2},..., \mathbf{R_M}, \mathbf{r_1}, \mathbf{r_2},..., \mathbf{r_N}, t)}{\partial t}= &\\
|
||
%
|
||
%&i\hbar\frac{\partial\Psi(\mathbf{R}_1,~\mathbf{R}_2,...,~ \mathbf{R}_M,~ \mathbf{r}_1,~ \mathbf{r}_2,...,~ \mathbf{r}_N,~ t)}{\partial t}= & \nonumber \\
|
||
%&\hat{H}(\mathbf{R}_1, \mathbf{R}_2,..., \mathbf{R}_M, \mathbf{r}_1, \mathbf{r}_2,..., \mathbf{r}_N, t)\Psi(\mathbf{R}_1, \mathbf{R}_2,..., \mathbf{R}_M, \mathbf{r}_1, \mathbf{r}_2,..., \mathbf{r}_N, t) &
|
||
%\end{eqnarray}
|
||
\begin{eqnarray}
|
||
\label{TDSE}
|
||
i\hbar\frac{\partial\Psi(\mathbf{R}_\alpha, \mathbf{r}_j, t)}{\partial t}=\hat{H}(\mathbf{R}_\alpha, \mathbf{r}_j, t)\Psi(\mathbf{R}_\alpha, \mathbf{r}_j, t), \quad \alpha=1, 2,..., M; j=1, 2,..., N.
|
||
\end{eqnarray}
|
||
This equation describes the evolution of the wavefunction in space and time, where $i$ is the imaginary unit. $\hbar$=$h/2\pi$ ~(1.054572$\times$ 10$^{-34}$~J$\cdot$s) is the reduced Planck constant and $t$ is the time. \textbf{R$_\alpha$} and \textbf{r}$_j$ refer to the coordinates of nucleus $\alpha$ and electron $j$, respectively. $\hat{H}$ is the so called Hamilton operator corresponding to the total energy of the system.
|
||
when the Hamiltonian itself is explicitly independent on time (the total wavefunction still has a time dependency), it is possible to decompose the space variables and the time variable to write the time-independent Schr{\"o}dinger equation:
|
||
\begin{eqnarray}
|
||
\label{TISE}
|
||
\hat{H}\Psi_k(\mathbf{R}_\alpha, \mathbf{r}_j)=E_k\Psi_k(\mathbf{R}_\alpha, \mathbf{r}_j)
|
||
\end{eqnarray}
|
||
where $E_k$ is the total energy of the system associated with the eigenstate $\Psi_k$. According to this, the evolution of the wavefunction becomes:
|
||
|
||
\begin{eqnarray}
|
||
\label{waveFunc}
|
||
\Psi(\mathbf{R}_\alpha, \mathbf{r}_j, t)=\sum_{k}c_k\Psi_k(\mathbf{R}_\alpha, \mathbf{r}_j)e^{-iE_kt/\hbar}
|
||
\end{eqnarray}
|
||
where $c_k$ is the coefficient at t=0. The eigenstates obtained from equation \ref{TISE} are the stationary states of the system and form a complete basis of
|
||
orthonormal vectors. The lowest energy eigenstate is called the ground state usually denoted as $\Psi_0$ and $E_0$ is the corresponding energy.
|
||
|
||
In a system made up of $M$ nuclei and $N$ electrons, the Hamiltonian operator (in a non-relativistic framework) is written as following:
|
||
%\begin{eqnarray}
|
||
\begin{align}
|
||
%\label{TISE}
|
||
\hat{H} &= \hat{T} + \hat{V} \nonumber \\
|
||
&=\hat{T}_n + \hat{T}_e + \hat{V}_{nn} + \hat{V}_{ee} + \hat{V}_{ne} \nonumber \\
|
||
&=-\frac{1}{2}
|
||
\sum_{\alpha=1}^{M}\frac{1}{M_\alpha}\mathbf{\nabla}_\alpha^2 - \frac{1}{2}
|
||
\sum_{j=1}^{N}\mathbf{\nabla}_j^2 + \sum_{\alpha=1}^{M}\sum_{\beta>\alpha}^{M}\frac{Z_\alpha Z_\beta}{|\mathbf{R}_\alpha - \mathbf{R}_\beta|} \nonumber \\
|
||
&\quad~ + \sum_{j=1}^{N}\sum_{i>j}^{N}\frac{1}{|\mathbf{r}_j - \mathbf{r}_i|} - \sum_{\alpha=1}^{M}\sum_{j=1}^{N}\frac{Z_\alpha}{|\mathbf{R}_\alpha-\mathbf{r}_j|}
|
||
\label{operatorH}
|
||
\end{align}
|
||
%\end{eqnarray}
|
||
in which the atomic unit system is used. The energy is thus expressed in Hartree. The mass is in unit of the mass of electron and the length is in Bohr. The vectorial variables are in bold in this present manuscript. $M_\alpha$ refers to the mass of nucleus $\alpha$ (in atomic unit) and $Z_\alpha$ is the atomic number of $a$. $\hat{T}_n$ and $\hat{T}_e$ are the kinetic energy operators of nuclei and electrons, respectively. $\hat{V}_{nn}$, $\hat{V}_{ee}$ and $\hat{V}_{ne}$ denote the potential energy operators of the repulsion between the nuclei, repulsion between the electrons and electrostatic attraction between nuclei and electrons, respectively.
|
||
\textbf{$\nabla^2$} is the Laplace operator. For nucleus $\alpha$ in Cartesian coordinates in three dimensions, its position vector is \textbf{R}$_\alpha$=(X$_\alpha$, Y$_\alpha$, Z$_\alpha$) and \textbf{$\nabla^2$} is expressed as following:
|
||
\begin{eqnarray}
|
||
\label{laplace}
|
||
\mathbf{\nabla}^2 = \frac{\partial^2}{\partial{X}_\alpha^2} +
|
||
\frac{\partial^2}{\partial{Y}_\alpha^2} + \frac{\partial^2}{\partial{Z}_\alpha^2}
|
||
\end{eqnarray}
|
||
|
||
\section{Born-Oppenheimer Approximation} \label{BO_approx}
|
||
|
||
Electrons are very light particles which can not be described correctly even qualitatively through classical mechanics. If we want to describe the electron
|
||
distribution in details, quantum mechanics must be applied, \textit{i.e.} solving the Schr{\"o}dinger equation. Usually in atomic and molecular systems,
|
||
it is very hard, if not impossible, to obtain the exact solutions of Schr{\"o}dinger equation, as this requires to separate variables. It is only possible in the
|
||
system containing one nucleus and one electron, \textit{i.e.}, the hydrogen atom or hydrogenic ions. For molecular species, the mathematical complexity
|
||
to solve the Schr{\"o}dinger equation will increase with the number of degrees of freedom of the system. Thus, it is necessary to resort to approximations
|
||
in almost all cases.
|
||
The \textbf{Born-Oppenheimer (BO) approximation} or adiabatic approximation is a cornerstone in real-life quantum analysis of atoms and molecules,
|
||
which helps to solve the Schr{\"o}dinger equation.
|
||
|
||
BO approximation is based on the large difference of mass between nuclei and electron and correspondingly the time scales of their motions.
|
||
Indeed, the mass of proton is much higher than that of electron ($m_p$/$m_e$ $\approx$ 1836). With the same amount of kinetic energy,
|
||
the electrons move much faster than the nuclei.
|
||
%($T_n$ \ll~ $T_e$).
|
||
%<<<<<<< HEAD
|
||
%Then, it is considered that the electrons move in the field of the fixed atomic nuclei, \textit{i.e.}, the electrons adapt instantly to the displacement of the nuclei while remaining in their ground state.\cite{Born1927}
|
||
|
||
%BO approximation consists of expressing the total wavefunction of a molecule as the product of a nuclear (vibrational, rotational) wavefunction and an electronic wavefunction, which enables a separation of Hamiltonian operator into the fast electronic term and the usually much slower nuclear term, where the coupling between electrons and nuclei is neglected so that the two smaller and decoupled systems can be solved more efficiently. In mathematical terms, the total wavefunction $\Psi_\mathrm{tot}$ of a molecule can be expressed as an expansion in the complete set of electronic wavefunctions $\psi_k^e$ with the expansion coefficients being functions of the nuclei coordinates $\psi^n$:
|
||
%=======
|
||
Then, it is considered that the electrons move in the field of the fixed atomic nuclei, \textit{i.e.}, the electrons adapt instantly to the displacement
|
||
of the nuclei while remaining in their ground state.\cite{Born1927}
|
||
|
||
BO approximation consists of expressing the total wavefunction of a molecule as the product of a nuclear (vibrational, rotational) wavefunction and
|
||
an electronic wavefunction, which enables a separation of the Hamiltonian operator into the fast electronic term and the usually much slower nuclear
|
||
term, where the coupling between electrons and nuclei is neglected so that the two smaller and non-coupled systems can be solved more efficiently.
|
||
In mathematical terms, the total wavefunction $\Psi_\mathrm{tot}$ of a molecule can be expressed as an expansion in the complete set of electronic
|
||
wavefunctions $\psi_k^e$ with the expansion coefficients being functions of the nuclei coordinates $\psi^n$:
|
||
\begin{eqnarray}
|
||
\label{BO}
|
||
\Psi_\mathrm{tot}(\mathbf{R}_\alpha, \mathbf{r}_j) = \sum_{k=1}^{\infty}
|
||
\psi_k^e (\mathbf{r}_j; \mathbf{R}_\alpha) \psi_k^n(\mathbf{R}_\alpha)
|
||
\end{eqnarray}
|
||
%<<<<<<< HEAD
|
||
%where the semicolon symbolizes the positions of the nuclei as parameters and not the variables of electronic wavefunction. This indicates although $\psi_k^e$ is a real-valued function of \textbf{r$_j$}, its functional form depends on \textbf{R$_\alpha$}.
|
||
%=======
|
||
where the semicolon symbolizes the positions of the nuclei as parameters and not the variables of electronic wavefunction.
|
||
This indicates that, although $\psi_k^e$ is a real-valued function of \textbf{r$_j$}, its functional form depends on \textbf{R$_\alpha$}.
|
||
|
||
Two smaller, consecutive steps can be used when using the BO approximation.
|
||
In the first step, the nuclei is treated as stationary.
|
||
%the nuclear kinetic energy is neglected,
|
||
The corresponding operator $\hat{T}_n$ is subtracted from the total molecular Hamiltonian operator $\hat{H}$ leading to the electronic Hamiltonian without considering nuclear kinetic energy:
|
||
\begin{eqnarray}
|
||
\label{H_e}
|
||
\hat{H}_e= \hat{T}_e + \hat{V}_{nn} + \hat{V}_{ee} + \hat{V}_{ne}
|
||
\end{eqnarray}
|
||
The major computational work is in solving the electronic Schr{\"o}dinger equation for a given set of nuclear coordinates.
|
||
\begin{eqnarray}
|
||
\label{BOSeq}
|
||
\hat{H}_{e}(\mathbf{R}_\alpha, \mathbf{r}_j)\psi_k^e(\mathbf{R}_\alpha, \mathbf{r}_j)=E_k^e\psi_k^e(\mathbf{R}_\alpha, \mathbf{r}_j)
|
||
\end{eqnarray}
|
||
where the eigenvalue $E_k^e$, electronic energy, depends on the chosen positions \textbf{R$_\alpha$} of the nuclei. By varying these positions \textbf{R$_a$} in small steps and repeatedly to solve the electronic Schr{\"o}dinger equation, so one can obtain $E_k^e$ as a function of \textbf{R$_\alpha$}, which generates the PES.
|
||
|
||
In the second step, the nuclear kinetic energy ${T}_n$ is reintroduced, and the Schr{\"o}dinger equation for the nuclear motion is:
|
||
\begin{eqnarray}
|
||
\label{nucSeq}
|
||
(\hat{T}_n + E_{k}^e(\mathbf{R}_\alpha) + \langle \Psi_k \vert \mathbf{\nabla}_n^2 \vert \Psi_k \rangle )\psi_k^n(\mathbf{R}_\alpha) = E_\mathrm{tot}\psi_k^n(\mathbf{R}_\alpha)
|
||
\end{eqnarray}
|
||
|
||
The eigenvalue $E_\mathrm{tot}$ is the total energy of the molecule, which includes the overall rotation translation of the molecule, contributions from electrons,
|
||
and nuclear vibrations. This second step involves a separation of vibrational, translational, and rotational motions.
|
||
Born and Oppenheimer assumed that the integral $\langle \Psi_k \vert \mathbf{\nabla}_n^2 \vert \Psi_k \rangle$ (diagonal correction) weakly depends on the nuclear coordinates,
|
||
so that it can be ignored. \cite{Born1927} Therefore, the Born-Oppenheimer approximation allows to describe the movement of nuclei in the corresponding potential to an adiabatic electronic state by the following equation:
|
||
\begin{eqnarray}
|
||
\label{ESeq}
|
||
(\hat{T}_n + E_{k}^e(\mathbf{R}_\alpha))\psi_k^n(\mathbf{R}_\alpha) = E_\mathrm{tot}\psi_k^n(\mathbf{R}_\alpha)
|
||
\end{eqnarray}
|
||
|
||
In this thesis, we will assume electrons adapt fast to reach their electronic ground state.
|
||
%I only focus on the electronic ground state.
|
||
The potential energy $E_0^e$ thus equals the ground state electronic
|
||
energy $E_{0}$ and the total energy is then $E_0^\mathrm{tot}=T_\mathrm{n} + E_{0}$. The Schr{\"o}dinger equation for $\psi_0^n(\mathbf{R}_\alpha)$ can
|
||
therefore be written as:
|
||
\begin{eqnarray}
|
||
\label{classical}
|
||
(\hat{T}_n + E_{0})\psi_0^n(\mathbf{R}_\alpha) = E_\mathrm{tot}\psi_0^n(\mathbf{R}_\alpha)
|
||
\end{eqnarray}
|
||
The next step is usually to consider the nuclei can be described classically. One can then consider that they evolve classically \textit{i.e.} following new Newton's equation, on a PES defined by the ground state electronic energy. The calculation of the ground state electronic energy is discussed in section~\ref{electronicEnergy}
|
||
%<<<<<<< HEAD
|
||
%The Born-Oppenheimer approximation only introduces very small errors for the majority of systems, which is widely applied in quantum chemistry to speed up the computation of molecular wavefunctions of large molecules. However, the BO approximation is only valid when the electronic state is sufficiently separated from other electronic states. That is to say when two or more solutions of the electronic Schr{\"o}dinger equation come close to each other energetically, the approximation loses validity (usually called ``break down") but it can be used as a starting point for more refined methods.\cite{Epstein1966, Butler1998}
|
||
%=======
|
||
The BO approximation only introduces very small errors in most systems, which explains it is widely applied in quantum chemistry
|
||
to speed up the computation of molecular wavefunctions of large molecules. However, the BO approximation is only valid when the electronic state
|
||
is sufficiently separated from other electronic states. That is to say that, when two or more solutions of the electronic Schr{\"o}dinger equation come
|
||
close to each other energetically, the approximation loses validity (usually called ``break down") but it can be used as a starting point for more refined
|
||
methods.\cite{Epstein1966, Butler1998}
|
||
%>>>>>>> 92023a10c3aa8b7dc4ace43987c1d571fb99a738
|
||
|
||
\section{Computation of Electronic Energy} \label{electronicEnergy}
|
||
|
||
%<<<<<<< HEAD
|
||
%The spectroscopic investigations in experiment assist to understand the electronic structure of molecules, for instance the measurements of absorption, radiation and scattering.
|
||
%These measurements often can provide a detailed picture of molecular systems but sometimes it is difficult to interpret.
|
||
%In the last few decades, molecular electronic-structure
|
||
%theory has developed to a stage where it can provide invaluable assistance in the interpretation of experimental measurements of a wide range important properties of molecules in rotational and vibrational spectroscopies, magnetic-resonance spectroscopies, and ultraviolet visible spectroscopies, etc.\cite{Puzzarini2010, Helgaker1999, Pedone2010, Fleming1986, Laane1999, Helgaker2012}
|
||
|
||
%The motion of systems including three or more masses interaction can not be solved analytically, so approximation methods must be applied. Many approximations have been proposed to obtain approximate solutions of the exact electronic wavefunction. Each of these approximations is usually the basis of one or more calculation approaches, which has their own advantages and disadvantages.
|
||
%When the generated solutions of the Schr{\"o}dinger equation without reference to experimental data, the methods are usually called \textit{ab initio} compared to semi-empirical models.
|
||
%Hartree-Fock (HF) model takes all interactions between electrons into account in an average fashion which means the average electron-electron interaction and the correlation between electrons are neglected.
|
||
%%In HF model, every electron is described by an orbital and the total wavefunction is given as a product of orbitals. HF equations depend on their own solutions which must be solved iteratively.
|
||
%Post HF theory usually generate more accurate calculation results with considering the electron correlation. \cite{Jensen2017}
|
||
%Semi-empirical methods are derived from HF model with neglecting all integrals involving more than two nuclei in the construction of Fock matrix to reduce the computational cost.
|
||
%Densituy functional theory (DFT) methods in the Kohn-Sham (KS) version can be regarded as an improvement on HF theory.
|
||
%Density Functional based Tight binding (DFTB) method is an abbreviation of DFT by doing additional approximations.\cite{dftb1, dftb2, Elstner1998, Elstner2014} Of course, there are premises and limitations of different methods.
|
||
%The PES can be explored with \textit{ab initio} methods. \cite{Cremer1993}
|
||
%Compared with quantum chemical (QM) methods that require considerable computer resources, molecular mechanics (MM) calculations are much cheaper but the types of reactions that can be dealed with MM method is very limited. The possibility of observing chemical reactions using \textit{ab initio} MD calculation is severely hindered by the short simulation time accessible.
|
||
%QM/MM techniques that combine QM for the reactive region and MM for the remainder are very promising, especially for big systems.\cite{Singh1986, Gao1996, Mordasini1998}
|
||
%\textit{ab initio} MD has been used to calculate the PES.\cite{Thompson1998, Schlegel2003}
|
||
|
||
|
||
%We will focus on the description of wavefunction based methods, density functional theory methods (DFT), and the density functional based tight-binding (DFTB) methods in this present section.
|
||
%a useful connection between density based and wavefunction based methods.
|
||
%\subsection{Wavefunction based Methods}
|
||
%\textbf{Hartree-Fock theory}
|
||
|
||
%HF theory is one of the earliest wavefunction based approximation methods, which is the foundation for a large part of the computational work on the electronic structure of atoms and molecules.\cite{Hartree1928, Hartree1947, Slater1968}
|
||
%The HF approximation transforms the many-body Schr{\"o}dinger equation into many coupled single-particle equations, which often assumes the $N$-body electronic wavefunction of a system is approximated by a single Slater determinant and every electron is considered to be independent.
|
||
%Hartree proposed that the electronic wavefunction could be approximated by assuming that in addition to the nuclei, the individual electrons could be separated as well. Therefore, the many-electron wavefunction would be a product of one-electron wavefunctions $\psi_j$:
|
||
%=======
|
||
Experimental spectroscopic investigations help in understanding the electronic structure of molecules, for instance measurements of absorption,
|
||
emission and scattering. These measurements can often provide a detailed picture of molecular systems but sometimes they are difficult to interpret.
|
||
In the last few decades, molecular electronic-structure theory has developed to a stage where it can provide invaluable
|
||
assistance in the interpretation of experimental measurements of a wide range of important properties of molecules in rotational and vibrational spectroscopies,
|
||
magnetic-resonance spectroscopies, ultraviolet/visible spectroscopies, and otherss.\cite{Puzzarini2010, Helgaker1999, Pedone2010, Fleming1986, Laane1999, Helgaker2012}
|
||
|
||
Electronic wavefunction of systems including three or more interacting particles can not be obtained analytically, so approximations must be applied.
|
||
Many approximations have been proposed to obtain approximate solutions of the exact electronic wavefunction. Each one of them is usually the basis
|
||
of one or more calculation approaches, which have their own advantages and disadvantages.
|
||
When solutions of the Schr{\"o}dinger equation are obtained without reference to experimental data, the methods are usually referred to as \textit{ab initio}
|
||
compared to semi-empirical models. The Hartree-Fock (HF) model takes all interactions between electrons into account except for the correlation between electrons
|
||
that is neglected.
|
||
Post-HF theory usually generates more accurate results by considering the electronic correlation.\cite{Jensen2017}
|
||
%Semi-empirical methods are derived from HF model with neglecting all integrals involving more than two nuclei in the construction of Fock matrix to reduce the computational cost.
|
||
DFT methods in the Kohn-Sham formulation can be regarded as an improvement over the HF theory as it considers
|
||
approximated electronic correlation.
|
||
Density-functional based tight-binding (DFTB) formalism is an approximated DFT method that involves additional approximations.\cite{dftb1, dftb2, Elstner1998, Elstner2014}
|
||
PES can be explored with \textit{ab initio} methods.\cite{Cremer1993}
|
||
Compared with quantum chemical (QM) methods that require considerable computer resources, molecular mechanics (MM) calculations are much cheaper but
|
||
present severe limitations in the treatment of chemical reactivity.. The possibility to model chemical reactions using \textit{ab initio} MD calculation is also
|
||
severely hindered by the short simulation time accessible with such methods. QM/MM techniques that combine QM for the reactive region and MM for the
|
||
remainder are very promising, especially for large systems.\cite{Singh1986, Gao1996, Mordasini1998}
|
||
|
||
The following section focuses on the description of wavefunction based methods, density functional theory method, and the density-functional based
|
||
tight-binding methods used to solve Schr{\"o}dinger equation and the electronic structure problem.
|
||
|
||
\subsection{Wavefunction based Methods}
|
||
|
||
\textbf{Hartree-Fock theory.}
|
||
Hartree-Fock theory is one of the earliest wavefunction based approximation methods, which is the foundation for a large part of the computational work on the electronic structure of atoms and molecules.\cite{Hartree1928, Hartree1947, Slater1968}
|
||
The HF approximation transforms the many-body Schr{\"o}dinger equation into many coupled single-particle equations, which often assumes the $N$-body electronic wavefunction of a system is approximated by a single Slater determinant and every electron is considered to be independent.
|
||
Hartree proposed that the electronic wavefunction could be approximated by assuming that in addition to the nuclei, the individual electrons could be separated as well. Therefore, the many-electron wavefunction would be a product of one-electron wavefunctions $\psi_j$:
|
||
|
||
%>>>>>>> 92023a10c3aa8b7dc4ace43987c1d571fb99a738
|
||
\begin{eqnarray}
|
||
\label{HF}
|
||
\Psi(\mathbf{r}_1, \mathbf{r}_2,..., \mathbf{r}_n)=\psi_1(\mathbf{r}_1)\psi_2(\mathbf{r}_2)...\psi_n(\mathbf{r}_n)
|
||
\end{eqnarray}
|
||
|
||
The HF theory assumes that every electron moves in an average field of all the other electrons and the nuclei in the molecule, which is an example of a mean-field approximation.
|
||
The HF equations for an individual electron $j$ moving in the mean field \textit{V}$_i^{HF}$, can be expressed:
|
||
|
||
\begin{equation}
|
||
\left(-\frac{1}{2}\mathbf{\nabla}_j^2 + V_\sigma(\mathbf{r}_j) + V_\mathrm{H}(\mathbf{r}_j)\right)\psi^\mathrm{HF}_\sigma (\mathbf{r}_j)-
|
||
\sum_{k=1}^{N_\sigma} \int d^3 {r'} \frac{\psi_{k\sigma}^{\mathrm{HF*}}(\mathbf {r'})\psi_{k\sigma}^\mathrm{HF}(\mathbf r_j)}{\mid \mathbf r_j-\mathbf{r'}\mid}\psi^\mathrm{HF}_{\sigma} (\mathbf{r'})= E_\sigma^\mathrm{HF} \psi_\sigma^\mathrm{HF}(\mathbf{r}_j)
|
||
\label{HF-individual}
|
||
\end{equation}
|
||
|
||
\begin{equation}
|
||
V_\mathrm{H}(\mathbf{r}_j) = \int d^3 {r}_j \frac{\rho (\mathbf{r'})}{\mid \mathbf{r}_j - \mathbf{r'} \mid}
|
||
\label{Hartreepotential}
|
||
\end{equation}
|
||
|
||
\begin{equation}
|
||
\rho(\mathbf{r}_j) =\sum_\sigma\sum_j^{N_\sigma}\rho_\sigma(\mathbf{r}_j)=\sum_\sigma\sum_j^{N_\sigma} \mid \psi_{\sigma}^\mathrm{H}(\mathbf {r}_j)\mid^2
|
||
\label{HF-idensity}
|
||
\end{equation}
|
||
%density means the probability of findind an electron
|
||
where $\nabla_j^2$ depends on the $j$th electron coordinates.
|
||
$V_\sigma(\textbf{r}_j)$ refers to the external potential. $\sigma$ is the spin. $V_{H}(\mathbf{r}_j)$ is the Hartree potential
|
||
% Hartree potential is a local potential.
|
||
%j\sigma refers to a state.
|
||
The last term on the left-hand side of eq \ref{HF-individual} is the HF exchange potential.
|
||
%which can be viewed as an operator acting on the ortital $\psi^{\mathbf{HF}_{j\sigma}(\mathbf r_j)$ in such a way that the result of its action at position $(\mathbf{r}_j}$ involves an integration of $\psi^{\mathbf{HF}_{j\sigma} over all \mathbf{r'}.
|
||
\textbf{r'} refers to the positions of all other electrons except electron $j$.
|
||
%The solutions of the individual Hartree equations depend on the solutions of all the other individual Hartree equations via \textit{V}$_i^{\mathrm{HF}}$. Thus, they have to be solved iteratively.
|
||
%<<<<<<< HEAD
|
||
%%First, an initial guess for the individual one-electron wavefunctions is needed. After this, iterate continues until all \psi$_{i}$ stay the same, which means a self-consistent field (SCF) has been achieved.
|
||
%The wavefunctions of the individual electrons are called orbitals. Because of the mean field, the movement of each electron is dependent on the orbits of the other electrons. The solutions of the single particle equation for atoms and molecules are the so called atomic orbitals (AOs) and molecular orbitals (MOs), respectively.
|
||
%HF or self-consistent field (SCF) theory can be derived by invoking the variational principle in a restricted space of wavefunctions, which was introduced in many books.\cite{Mcweeny1992, Szabo1996, Kohanoff2006, Helgaker2014}
|
||
%Molecular orbitals are polycentric, which is difficult to solve the HF equations of molecules.
|
||
%Roothaan put forward the approximation of MOs as a linear combination of atomic orbitals (LCAO) in 1951, namely, a linear combination of atomic basis functions to solve the HF equations of molecules, which is a significant improvement in the practical solutions of the HF equations. \cite{Roothaan1951} Roothaan equations transform a HF computation to a problem in matrix algebra, for which the algebraic equations are particularly suitable for modern computers. \cite{Blinder2018} Both the vast majority of the calculation approaches, for instance, $ab \ initio$, density functional, semi-empirical molecular orbital, or even some sophisticated force-fields and the qualitative comprehension of chemistry are based on the concept that the orbitals of a given molecule can be built from the orbitals of the constituent atoms. \cite{koch1999}
|
||
|
||
%The HF theory showed some success. \cite{Moskowitz1965}
|
||
%However, it does not contain the electron correlation beyond the minimum required to satisfy the antisymmetry for electronic wavefunctions, and the resulting approximate electronic energies are not accurate enough for most practical applications in chemistry.
|
||
%Good HF results can account for over 99\% of the true total energy of the system. Due to the variational principle, the HF wavefunction is too high in energy. So the remaining 1\% error of the true total energy which is defined as correlation energy ($E_\mathrm{corr} = E_\mathrm{exact} - E_\mathrm{HF}$) is essential to account for the chemical properties of atoms and molecules. Electron correlation is from the correlated movement of electrons, and the failure of HF theory to describe it originates from that the mean-field approximation can not treat electron-electron interactions properly.\cite{Lowdin1995}
|
||
|
||
%\textbf{Post Hartree-Fock theory}
|
||
|
||
%Post-Hartree-Fock methods are the improvements of HF, which add the electron correlation.\cite{Cramer2003, Jensen2017} The calculation of the correlation energy is then the objective of several post-HF methods (for instance, the configuration interaction (CI) method, \cite{Head1994, Maurice1999} Møller-Plesset perturbation theory (MP2, MP3 and MP4) methods, \cite{Moller1934, Pople1976, Krishnan1978} coupled cluster (CC) method, \cite{Vcivzek1966, Purvis1982, Raghavachari1989, Van2001} quantum chemistry composite (G2, G3, and T1) methods \cite{Curtiss1991, Curtiss1998, Ohlinger2009} and so on). The post-Hartree-Fock methods \cite{Head1994, Magnasco2009, Dacosta2011} usually give more accurate results than the HK calculations, but the added accuracy comes with the price of more computational cost. Among the above post-HF methods, the MP2 method was used in this thesis.
|
||
|
||
%The major calculation of HF is the variational optimization of the orbitals of a single configuration state function (CSF). In the case of a closed-shell singlet and some open-shell cases, this would be a single Slater determinant. To include the electron correlation in wavefunction theory, a superposition of CSFs can be considered. A multiconfigurational wavefunction is a linear combination of two or more CSFs. Standard multiconfigurational approaches are the CI, CC, and perturbation theory methods.
|
||
|
||
|
||
%\subsection{Density Functional Theory}
|
||
%For a long time, approximations based on the wavefunction theory were typically applied to solve the Schr{\"o}dinger equation. However, it is usually impractical to perform a wavefunction theory calculation with chemical accuracy for complex or large systems. Density functional theory is based on electron density rather than electronic wavefunctions.\cite{Kohn1965, kohn1999nobel}
|
||
%Because of having more favorable scaling of demand on computational resources with respect to system size, DFT is the most widely used method available in computational chemistry, computational physics, and condensed-matter physics for ground state of large and complex systems.
|
||
%DFT makes it possible to transform the problem of electrons interacting evolving in a nuclear potential to a problem of independent electrons in an effective potential.
|
||
%The electron density $\rho$(\textbf{r}) corresponds to the number of electrons per unit volume in a given state.
|
||
%The central idea of DFT is to promote the $\rho$($\mathbf{r}$) (functional which only depends on spatial coordinates) as a key variable in the determination of the energy of a system.
|
||
%From the electronic wavefunction of the system, $\rho$($\mathbf{r}$) is written as:
|
||
%=======
|
||
|
||
%First, an initial guess for the individual one-electron wavefunctions is needed. After this, iterate continues until all \psi$_{i}$ stay the same, which means a self-consistent field (SCF) has been achieved.
|
||
%The wavefunctions of the individual electrons are called orbitals. Because of the mean field, the movement of each electron is dependent on the orbits of the other electrons. The solutions of the single particle equation for atoms and molecules are the so called atomic orbitals (AOs) and molecular orbitals (MOs), respectively.
|
||
HF or self-consistent field (SCF) theory can be derived by invoking the variational principle in a restricted space of wavefunctions,
|
||
which was introduced in many books.\cite{Mcweeny1992, Szabo1996, Kohanoff2006, Helgaker2014}
|
||
C. Roothaan then put forward the approximation of molecular orbitals (MOs) as a linear combination of atomic orbitals (LCAO) in 1951, namely,
|
||
a linear combination of atomic basis functions to solve the HF equations of molecules.\cite{Roothaan1951} Thus was a significant improvement
|
||
in the practical solution of the HF equations. Roothaan equations allow to transform the HF problem into a linear algebra problem, for which
|
||
algebraic equations are particularly suitable for modern computers.\cite{Blinder2018} The vast majority of computational approaches, whether $ab \ initio$,
|
||
semi-empirical, or even some sophisticated force-fields are based on the concept that the molecular orbitals of a given molecule can be built
|
||
from the atomic orbitals of its constituent atoms.\cite{koch1999}
|
||
|
||
The HF theory showed some success.\cite{Moskowitz1965} However, it does not contain the electron correlation beyond the minimum required
|
||
to satisfy the antisymmetry for electronic wavefunctions. The resulting approximated electronic energies are therefore not accurate enough for
|
||
most practical applications in chemistry. Good HF results can account for over 99\% of the true total energy of the system. Due to the variational
|
||
principle, the HF wavefunction is always too high in energy. So, the remaining 1\% error with respect to the true total energy, which is defined as
|
||
the correlation energy ($E_\mathrm{corr} = E_\mathrm{exact} - E_\mathrm{HF}$) is essential to account for the chemical properties of atoms and
|
||
molecules. Electron correlation results from the correlated behavior of electrons, and the failure of HF theory to describe it originates from that
|
||
the mean-field approximation can not treat electron-electron interactions properly.\cite{Lowdin1995}
|
||
|
||
\textbf{Post Hartree-Fock methods.}
|
||
Post Hartree-Fock methods provide improvements to HF theory by adding the electron correlation.\cite{Cramer2003, Jensen2017}
|
||
The calculation of the correlation energy is then the objective of several post-HF methods (for instance, the configuration interaction
|
||
(CI) method,\cite{Head1994, Maurice1999} Møller-Plesset perturbation theory (MP2, MP3 and MP4) methods,\cite{Moller1934, Pople1976, Krishnan1978}
|
||
coupled cluster (CC) method,\cite{Vcivzek1966, Purvis1982, Raghavachari1989, Van2001} quantum chemistry composite (G2, G3, and T1) methods,\cite{Curtiss1991, Curtiss1998, Ohlinger2009} and so on). Post-Hartree-Fock methods usually give more accurate results than Hartree-Fock calculations, \cite{Head1994, Magnasco2009, Dacosta2011}
|
||
but the additional accuracy comes to the price of a higher computational cost. Among the aforementioned post-HF methods,
|
||
the Møller-Plesset perturbation theory at second-order (MP2) method was used along this thesis.
|
||
|
||
HF theory variationally optimize the orbitals of a single configuration state function (CSF). In the case of a closed-shell singlet and some
|
||
open-shell cases, this would be a single Slater determinant. To include the electron correlation in wavefunction theory, a superposition
|
||
of CSFs should be considered. This superposition, referred to as a multiconfigurational wavefunction, is a linear combination of two or
|
||
more CSFs. Standard multiconfigurational approaches are the CI, CC, and perturbation theory methods.
|
||
|
||
\subsection{Density Functional Theory}
|
||
For a long time, approximations based on the wavefunction were systematically applied to solve the Schr{\"o}dinger equation.
|
||
However, it is usually impractical to perform a wavefunction based calculation with chemical accuracy for complex or large systems.
|
||
Density functional theory is based on the electron density rather than the electronic wavefunction.\cite{Kohn1965, kohn1999nobel}
|
||
Because DFT displays a more favourable scaling of computational resources with respect to system size, DFT is nowadays the most
|
||
widely used method available in computational chemistry, computational physics, and condensed-matter physics for ground state
|
||
calculation of large and complex systems.
|
||
|
||
DFT makes it possible to transform the problem of electrons interacting and evolving in a nuclear potential to a problem of
|
||
independent electrons evolving in an effective potential. The electron density $\rho$(\textbf{r}) corresponds to the number of
|
||
electrons per unit volume in a given state.
|
||
The central idea of DFT is to promote $\rho$($\mathbf{r}$) (function which only depends on three spatial coordinates) as the
|
||
key variable in the determination of the electronic energy of a system. From the electronic wavefunction of the system, $\rho$($\mathbf{r}$)
|
||
is written as:
|
||
|
||
%>>>>>>> 92023a10c3aa8b7dc4ace43987c1d571fb99a738
|
||
\begin{eqnarray}
|
||
\label{density}
|
||
\rho(\mathbf{r})=N\int \Psi^*(\mathbf{r}, \mathbf{r}_2, \mathbf{r}_3, ..., \mathbf{r}_N) \Psi(\mathbf{r}, \mathbf{r}_2, \mathbf{r}_3,..., \mathbf{r}_N) d\mathbf{r}_2 d\mathbf{r}_3...d\mathbf{r}_N
|
||
\end{eqnarray}
|
||
|
||
This idea originates from the the model of the uniform electron gas in the phase space around an atom developed in 1927 by Thomas \cite{Thomas1927} and Fermi \cite{Fermi1928}, which is the predecessor to density functional theory. Nevertheless, the Thomas-Fermi model is unable to correctly describe molecular bonds because it does not take into account the exchange and correlation energies.
|
||
|
||
Although DFT has a history almost as old as the Schr{\"o}dinger equation, the modern form dates back to the paper published by P. Hohenberg and W. Kohn \cite{Hohenberg1964}
|
||
that introduced the two Hohenberg–Kohn (HK) theorems in 1964 and the extension by Levy in 1979.\cite{Levy1979} The theory is usually applied in the form latter suggested
|
||
by Kohn and Sham in 1965.\cite{Kohn1965}
|
||
|
||
The first HK theorem shows that, for a many-electron system in its ground state, the energy is uniquely determined by the electron density $\rho$(\textbf{r}).
|
||
In other words, the first HK theorem shows that it is not necessary to know the wavefunction of the system to obtain its energy and that the knowledge of the electron
|
||
density alone is sufficient. It sets down the foundation for reducing the many-body problem of $N$ electrons with $3N$ spatial coordinates to three spatial
|
||
coordinates, through using functionals of the electron density. From this theorem, it follows that $\rho$(\textbf{r}) determines the external potential $V_\mathrm{ext}$(\textbf{r})
|
||
and $N$ can be obtained via the normalization of $\rho$(\textbf{r}):
|
||
|
||
%<<<<<<< HEAD
|
||
%The first HK theorem shows that, for a many-electron system in its ground state, the energy is uniquely determined by the electron density $\rho$(\textbf{r}). In other words, HK theorem shows that it is not necessary to know the wavefunctional of the system to obtain its energy and that the knowledge of the electron density alone is sufficient. It set down the foundation for reducing the many-body problem of $N$ electrons with $3N$ spatial coordinates to three spatial coordinates, through using functionals of the electron density.
|
||
%From the first Hohenberg-Kohn (HK) theorem,\cite{Hohenberg1964} the $\rho$(\textbf{r}) determines the external potential $V_\mathrm{ext}$(\textbf{r}) and $N$ can be obtained via the normalization of $\rho$(\textbf{r}):
|
||
%=======
|
||
%>>>>>>> 92023a10c3aa8b7dc4ace43987c1d571fb99a738
|
||
\begin{eqnarray}
|
||
\label{N}
|
||
\int \rho(\mathbf{r}) d\mathbf{r}=N
|
||
\end{eqnarray}
|
||
and $N$ and $V_\mathrm{ext}$(\textbf{r}) determine the electronic Hamiltonian. $\rho$(\textbf{r}) determines the energy and all other ground state electronic properties of a system. This is clearly shown in Figure \ref{variables}.
|
||
\figuremacro{variables}{Interdependence of basic variables in the Hohenberg-Kohn theorem.} \\
|
||
|
||
The second HK theorem is a variational electron density theorem which defines an energy functional for a system.
|
||
For a given external potential $V_\mathrm{ext}$(\textbf{r}), the ground state energy $E_0$ of the system is the minimum
|
||
value of the energy functional obtained for the exact ground state electron density $\rho$(\textbf{r}) of the system.
|
||
%%
|
||
From the HK theorems, we can write the functional of the total energy of the system as a sum of the kinetic energy of th
|
||
electrons $T_\mathrm{e}$[$\rho$(\textbf{r})], and the electronic interaction energy $E_\mathrm{ee}$[$\rho$(\textbf{r})]:
|
||
\begin{eqnarray}
|
||
\label{HK}
|
||
E[\rho(\mathbf{r})]=T_\mathrm{e}[\rho(\mathbf{r})]+V_\mathrm{ee}[\rho(\mathbf{r})]+ \int V_\mathrm{ext}(\textbf{r})\rho(\mathbf{r})d(\mathbf{r})
|
||
\end{eqnarray}
|
||
the last term represents the interaction between the lectron density and the external potential, \textit{i.e.} the nuclei in the case of atoms and molecules.
|
||
To obtain the energy of the ground state from this equation, the variation principle can be applied with respect to $\rho$(\textbf{r}). To do this, the form of $T_\mathrm{e}$[$\rho$(\textbf{r})] and $V_\mathrm{ee}$[$\rho$(\textbf{r})] should be known.
|
||
|
||
The HK theorems do not provide mathematical expressions for $T_\mathrm{e}$[$\rho$(\textbf{r})] and $V_\mathrm{ee}$[$\rho$] for a system
|
||
of interacting particles. To solve this problem, Kohn and Sham first proposed to treat the electrons as non-interacting particles subject to the
|
||
$V_\mathrm{ext}$(\textbf{r}) potential only.\cite{Kohn1965} The idea was to work with a fictitious system of $N$ non-interacting
|
||
electrons evolving in an effective potential and having the exact electronic density of the system. The resolution of equation~\ref{HK}
|
||
for a system of non-interacting electrons is known exactly, and if the correct electron density is reproduced, then the exact electronic
|
||
energy of the system can be calculated. The total energy of the real system is composed as follows:
|
||
|
||
\begin{eqnarray}
|
||
\label{E_KS}
|
||
E_\mathrm{DFT}[\rho(\mathbf{r})]=T_\mathrm{no}[\rho(\mathbf{r})]+\underbrace{E_\mathrm H[\rho(\textbf{r})] + \int V_\mathrm{ext}(\textbf{r})\rho(\mathbf{r})d(\mathbf{r}) +E_\mathrm{xc}[\rho(\mathbf{r})]}_{E_\mathrm{eff}[\rho(\mathbf{r})]}
|
||
\end{eqnarray}
|
||
where the sum of last three terms on the right-hand side is the effective energy $E_\mathrm{eff}$[$\rho(\textbf{r}$)]. \\
|
||
$T_\mathrm{no}$[$\rho$(\textbf{r})] is the kinetic energy of a system of non-interacting electrons:
|
||
|
||
\begin{eqnarray}
|
||
\label{Tn}
|
||
T_\mathrm{no}[\rho(\mathbf{r})]=\sum_{i}^{N} \left \langle \Psi_i \left|-\frac{1}{2}\nabla^2 \right | \Psi_i \right \rangle
|
||
\end{eqnarray}
|
||
|
||
$E_H$[$\rho$(\textbf{r})] represents the Hartree energy which corresponds to the interaction energy of a classical charge distribution
|
||
of density $\rho$(\textbf{r}):
|
||
|
||
\begin{eqnarray}
|
||
\label{Ehartree}
|
||
E_\mathrm{H}[\rho(\textbf{r})]=\frac{1}{2} \int \int \frac{\rho(\mathbf{r})\rho(\mathbf{r'})}{|\mathbf{r}-\mathbf {r'}|}d\mathbf{r}d\mathbf{r'}
|
||
\end{eqnarray}
|
||
|
||
$V_\mathrm{ext}$(\textbf{r}) is the external potential.
|
||
The remaining energy components are assembled in the exchange-correlation energy $E_\mathrm{xc}$[$\rho$(\textbf{r})] functional containing
|
||
the difference between the kinetic energy of the real system $T$[$\rho$(\textbf{r})] and that of the non-interacting system $T_n$[$\rho$(\textbf{r})]
|
||
and the non-classical part of $E_\mathrm{ee}$[$\rho$].
|
||
The $E_\mathrm{xc}$[$\rho$(\textbf{r})] functional can thus be expressed as:
|
||
\begin{eqnarray}
|
||
\label{Exc}
|
||
E_\mathrm{xc}[\rho(\mathbf{r})]=(T[\rho(\mathbf{r})]-T_\mathrm{no}[\rho(\mathbf{r})])+(V_\mathrm{ee}[\rho(\mathbf{r})]-E_\mathrm H[\rho(\textbf{r})])
|
||
\end{eqnarray}
|
||
|
||
To minimize the energy $E$[$\rho$(\textbf{r})] with respect to $\rho$(\textbf{r}) applying by the variational principle
|
||
while considering the constraints of orbital orthogonality, one performs an optimization under the constraints using Lagrange multipliers.
|
||
For the optimal $\rho$(\textbf{r}), the energy $E$[$\rho$(\textbf{r})] does not change upon the variation of $\rho$($\mathbf{r}$), providing that
|
||
$\rho$(\textbf{r}) obey to equation~\ref{N}:
|
||
|
||
\begin{eqnarray}
|
||
\label{Delta}
|
||
\delta(E[\rho(\mathbf{r})]-\mathcal{L}\rho(\mathbf{r}))=0
|
||
\end{eqnarray}
|
||
where $\mathcal{L}$ is a corresponding Lagrange multiplier.
|
||
Combining eqs \ref{E_KS}, \ref{Tn}, \ref{Ehartree}, and \ref{Delta}, the Euler equation can be written as:
|
||
|
||
\begin{eqnarray}
|
||
\label{Lagrange}
|
||
\mathcal{L}=V_\mathrm{eff}[\rho(\mathbf{r})]+ \frac{\delta T_n[\rho(\mathbf{r})]}{\delta \rho(\mathbf{r})}
|
||
\end{eqnarray}
|
||
where the effective potential $V_\mathrm{eff}$[$\rho$($\mathbf{r}$)] is introduced:
|
||
|
||
\begin{align}
|
||
V_\mathrm{eff}[\rho(\mathbf{r})] &=V_\mathrm{ext}[\rho(\mathbf{r})]+V_\mathrm{H}[\rho(\mathbf{r})]+V_\mathrm{xc}[\rho(\mathbf{r})] \nonumber \\
|
||
&=V_\mathrm{ext}[\rho(\mathbf{r})]+\int\frac{\rho(\mathbf{r'})}{|\mathbf{r}-\mathbf {r'}|}d\mathbf{r'} +\frac{\delta E_\mathrm{xc}[\rho(\mathbf{r})]}{\delta \rho(\mathbf{r})}
|
||
\label{Veff}
|
||
\end{align}
|
||
where $V_{H}$[$\rho$(\textbf{r})] refers to the Hartree potential and $V_{xc}$[$\rho$(\textbf{r})] is the exchange-correlation potential.
|
||
Therefore, for a given $V_\mathrm{eff}$[$\rho$(\textbf{r})], one can obtain $\rho$(\textbf{r}) by making the right-hand side of equation~\ref{Lagrange}
|
||
independent of \textbf{r} which is done by introducing the molecular orbitals $\psi_{i}$(\textbf{r}) such that:
|
||
|
||
\begin{eqnarray}
|
||
\label{rho}
|
||
\rho(\mathbf{r})=\sum_{i}^{N} \left |\psi_i(\mathbf{r}) \right |^2
|
||
\end{eqnarray}
|
||
%<<<<<<< HEAD
|
||
%where \textbf{x} denotes the four vector-containing space and spin variables, and the integration is implemented over the spin variable $\sigma$. Moreover, the molecular orbitals $\psi_{i}$[$\rho$(\textbf{r})] should satisfy the one-electron KS equation,
|
||
%\begin{eqnarray}
|
||
%\label{KS}
|
||
%\underbrace{\left(-\frac{1}{2}\mathbf{\nabla}_i^2+V_\mathrm{eff}[\rho(\mathbf{r})]\right)}_\mathrm{KS ~ operator} \psi_i[\rho(\mathbf{r})]=E_i\psi_i[\rho(\mathbf{r})]
|
||
%=======
|
||
|
||
The molecular orbitals $\psi_{i}$(\textbf{r}) should satisfy the one-electron KS equations:
|
||
|
||
\begin{eqnarray}
|
||
\label{KS}
|
||
\underbrace{\left(-\frac{1}{2}\mathbf{\nabla}_i^2+V_\mathrm{eff}[\rho(\mathbf{r})]\right)}_\mathrm{KS ~ operator} \psi_{i}(\textbf{r})=E_i \psi_{i}(\textbf{r})
|
||
%>>>>>>> 92023a10c3aa8b7dc4ace43987c1d571fb99a738
|
||
\end{eqnarray}
|
||
|
||
This result can be reobtained within a variational context when looking for those orbitals minimizing the
|
||
energy functional of equation~\ref{HK}, subject to orthonormality conditions:
|
||
\begin{eqnarray}
|
||
\label{orthonormality}
|
||
\int \psi^*_i(\mathbf{r}) \psi_j(\mathbf{r}) d\mathbf{r}=\delta_{ij}
|
||
\end{eqnarray}
|
||
|
||
%<<<<<<< HEAD
|
||
%The one-electron KS equation (eq \ref{KS}), just as the Hartree and HF equations, needs to be solved iteratively.
|
||
%The cost for the incorporation of electron correlation is to calculate the $V_{xc}$[$\rho$(\textbf{r})].
|
||
%=======
|
||
The one-electron KS equation of equation~\ref{KS}, just as the HF equations, needs to be solved iteratively.
|
||
The computational cost for the incorporation of electron correlation is the one necessary to calculate $V_{xc}$[$\rho$(\textbf{r})].
|
||
%>>>>>>> 92023a10c3aa8b7dc4ace43987c1d571fb99a738
|
||
%%the appearance of the exchange correlation potential, $V_\mathrm{xc}$[$\rho$(\textbf{r})].
|
||
|
||
DFT is in principle correct if one knows the “exact” exchange-correlation functional. However, despite a lot of work aimed at determining
|
||
this exact functional, its form is still unknown and a systematic strategy for improvement is not available. Therefore, it is always necessar
|
||
to use an approximate functional that is characterized by more or less important artefacts in effective calculations. Over the last decades,
|
||
many exchange-correlation functionals have been proposed. Although they are different, it is possible to classify them into families according
|
||
to some common characters. \textbf{Local Density Approximation} (LDA), \textbf{Generalized Gradient Approximation} (GGA), \textbf{meta-GGA},
|
||
and \textbf{hydrid functionals} (containing to some extent a contribution of HF exact-exchange) are some of the most widely used approximations.
|
||
\textbf{LDA} is the first approximation of $E_{xc}$[$\rho$(\textbf{r})] proposed by Kohn and Sham,\cite{Kohn1965} which is based on the description
|
||
of the homogeneous electron gas. For atomic or molecular systems, when densities vary rapidly in space, the assumption of a uniform electron
|
||
density is not correct and the LDA approximation is not applicable any more. For example, binding energies in molecules are usually overestimated
|
||
by the LDA approximation. To solve this problem, a new family of functionals, \textbf{GGA functional}, which includes a contribution of the electron
|
||
density gradient was developed.\cite{Kohn1965, Herman1969, Perdew1996}
|
||
\textbf{Hybrid functionals} were first proposed by Becke in 1993.\cite{Becke1993} The main idea is that for an uncorrelated system, the HF energy
|
||
is exact while for a highly correlated system LDA or GGA energy is more appropriate. These two states are then connected by a continuum of
|
||
partially correlated real systems of identical density. This connection is described by the adiabatic connection formula:
|
||
|
||
\begin{eqnarray}
|
||
\label{Ehybrid}
|
||
E_\mathrm{xc}[\rho(\mathbf{r})]=\int_{0}^{1} E_\mathrm{xc}[\rho(\mathbf{r})] \lambda d\lambda
|
||
\end{eqnarray}
|
||
%<<<<<<< HEAD
|
||
%where $\lambda$=0 corresponds to the uncorrelated case and $\lambda$=1 corresponds to the highly correlated case. Hybrid functionals usually provide higher quality results than those provided by LDA and GGA functionals for the study of molecular properties.
|
||
%Hybrid functionals have been very successful for ground-state properties. B3LYP (Becke, 3-parameter, Lee-Yang-Parr),\cite{Becke1988, Lee1988, Becke1993, Stephens1994} PBE with one empirical parameter (PBE0, also called PBE1PBE),\cite{Adamo1999} Heyd-Scuseria-Ernzerhof (HSE06),\cite{Henderson2009, Krukau2006}, and M06-2X \cite{Zhao2008} are some of the most popularly used Hybrid functionals.
|
||
%Approximate DFT suffers from some generic problems, many of which can be traced to the so-called delocalization error \cite{Li2018} or the closely related problem of self-interaction error.
|
||
%Nevertheless, it should be emphasized that today DFT, cast in the KS formalism, provides a computational tool with an remarkable quality and computationally less expensive than wavefunction theory based methods.
|
||
%=======
|
||
where $\lambda$=0 corresponds to the uncorrelated case and $\lambda$=1 corresponds to the highly correlated case.
|
||
Hybrid functionals usually provide higher quality results than those provided by LDA and GGA functionals for the study of
|
||
molecular properties.
|
||
Hybrid functionals have been very successful for ground-state properties. B3LYP (Becke, 3-parameter, Lee-Yang-Parr),\cite{Becke1988, Lee1988, Becke1993, Stephens1994}
|
||
PBE with one empirical parameter (PBE0, also called PBE1PBE),\cite{Adamo1999} Heyd-Scuseria-Ernzerhof (HSE06),\cite{Henderson2009, Krukau2006}, and
|
||
M06-2X \cite{Zhao2008} are some of the most popularly hybrid functionals.
|
||
Approximate functionals suffer from some generic problems, many of which can be traced to the so-called delocalization error \cite{Li2018} or
|
||
the closely related problem of self-interaction error. Nevertheless, it should be emphasized that today DFT, cast in the Kohn-Sham formalism,
|
||
provides a computational tool with a remarkable quality and computationally less expensive than wavefunction based methods.
|
||
%>>>>>>> 92023a10c3aa8b7dc4ace43987c1d571fb99a738
|
||
|
||
\subsection{Density Functional based Tight-Binding Theory} \label{sec:DFTB}
|
||
|
||
DFT is computationally too expensive for systems with more than hundreds of atoms especially when one needs to perform global optimization
|
||
or to perform molecular dynamics (MD) simulations with sufficient statistical sampling of the initial conditions or to perform fairly long trajectories.
|
||
It is therefore necessary to further simplify the method in order to reduce the computational cost. To do so, Foulkes and Haydock showed that
|
||
tight-binding models can be derived from the DFT.\cite{Foulkes1989} Later, DFTB was proposed by D. Porezag $et ~ al.$.\cite{Porezag1995}
|
||
|
||
Non-self-consistent DFTB scheme is suitable to study systems in which polyatomic electronic density is well described by a sum of atom-like densities.
|
||
This is the case ighly ionic and homonuclear covalent systems. However, uncertainty rises in the non-self-consistent DFTB scheme whem chemical bonds
|
||
are controlled by a subtler charge balance between atoms especially for polar, semi-conductor and heteronuclear molecules. The self-consistent-charge
|
||
extension of DFTB, SCC-DFTB, was born as an improvement of standard DFTB to provide a better description of electronic systems in which long-range
|
||
Coulomb interactions are significant.\cite{Elstner1998, Porezag1995, Seifert1996, Elstner1998, Elstne2007, Oliveira2009} However, this method is not
|
||
fully suitable for the calculation of molecular aggregates because of the dispersion and charge resonance present in such systems. Further corrections
|
||
were made to the model in order to address these problems. I know describe in more details how SCC-DFTB can be derived from DFT.
|
||
|
||
\textbf{General principles of DFTB.} DFTB is derived from DFT based on the following approximations:
|
||
|
||
\begin{itemize}
|
||
|
||
\item[$\bullet$] Only valence electrons are treated explicitly.
|
||
|
||
\item[$\bullet$] Molecular orbitals are developed on atomic valence orbitals.
|
||
|
||
\item[$\bullet$] A Taylor expansion of the total energy around a reference density is realized.
|
||
|
||
\item[$\bullet$] Integrals involving more than two centers are neglected.
|
||
\end{itemize}
|
||
|
||
%<<<<<<< HEAD
|
||
%The first order DFTB (historically called the zeroth order DFTB) takes into account the first term of Taylor's expansion, which is an equivalence with the method of tight-binding. The second order DFTB (historically the SCC-DFTB) is a self-consistent calculation on atomic charges. There is also a recent third order extension of the DFTB (called the DFTB3) \cite{Gaus2011} which is not used in this thesis so the third order extension part will not be shown in the following equations. According to eq \ref{E_KS}, $E_\mathrm{DFT}[\rho(\mathbf{r})]$ can be written as:
|
||
%=======
|
||
\textbf{First-order DFTB} (historically referred to as zeroth-order DFTB) takes into account the first term of Taylor's expansion and which is
|
||
equivalence with the other tight-binding model. \textbf{Second-order DFTB} (historically the SCC-DFTB) introduces a self-consistent procedure on
|
||
atomic charges. There is also a more recent \textbf{third-order extension of DFTB} (referred to as DFTB3).\cite{Gaus2011} DFTB3 was not used in
|
||
this thesis so the third-order expansion term will not be shown in the following equations. According to equation~\ref{E_KS}, $E_\mathrm{DFT}[\rho(\mathbf{r})]$
|
||
can be written as:
|
||
%>>>>>>> 92023a10c3aa8b7dc4ace43987c1d571fb99a738
|
||
\begin{align}
|
||
E_\mathrm{DFT}[\rho(\mathbf{r})]=&\langle \psi_i \vert \hat{T}_\mathrm{no}[\rho(\mathbf{r})]\vert \psi_i \rangle + E_\mathrm{eff}[\rho(\mathbf{r})]
|
||
\label{Edft}
|
||
\end{align}
|
||
|
||
The SCC-DFTB scheme is built on the second-order Taylor series development around a reference electronic density $\rho_{_0}$(\textbf{r}), so
|
||
$\rho$(\textbf{r})=$\rho_{_0}$(\textbf{r})+$\delta \rho$(\textbf{r}). In practice, $\rho_{_0}$(\textbf{r}) is taken as the superimposition of the
|
||
densities of isolated atoms. Then $E_\mathrm{DFTB}[\rho_{_0}(\mathbf{r})+\delta \rho(\mathbf{r})]$ can be written as follows:
|
||
|
||
\begin{align}
|
||
E_\mathrm{DFTB}[\rho_{_0}(\mathbf{r})+\delta \rho(\mathbf{r})]
|
||
=&\langle \psi_i \vert \hat{T}_\mathrm{no}[\rho(\mathbf{r})]\vert \psi_i \rangle + E_{eff}[\rho_{_0}(\mathbf{r})]+\int \frac{\delta E_{eff}[\rho(\mathbf{r})]}{\delta\rho(\mathbf{r})}\Big\vert _{\rho_{_0}}\delta \rho(\mathbf{r})d\mathbf{r} \nonumber \\
|
||
& +\frac{1}{2} \iint\frac{\delta^2 E_{eff}[\rho(\mathbf{r})]}{\delta\rho(\mathbf{r}) \delta\rho(\mathbf{r'})}\Big\vert _{\rho_{_0}, ~ \rho_{_0}'} \delta\rho(\mathbf{r}) \delta\rho(\mathbf{r'})d\mathbf{r}d\mathbf{r'}
|
||
\label{Edftb}
|
||
\end{align}
|
||
with $\delta \rho$(\textbf{r}) = $\rho$(\textbf{r})-$\rho_{_0}$(\textbf{r}). Equation~\ref{Edftb} can be rewritten as:
|
||
|
||
\begin{align}
|
||
E_\mathrm{DFTB}[\rho_{_0}(\mathbf{r})+\delta \rho(\mathbf{r})]=& \overbrace{E_\mathrm{eff}[\rho_{_0}(\mathbf{r})]-\int \frac{\delta E_\mathrm{eff}[\rho(\mathbf{r})]}{\delta\rho(\mathbf{r})}\Big\vert _{\rho_{_0}} \rho_{_0}(\mathbf{r})d\mathbf{r}}^{E_\mathrm{rep}} \nonumber \\
|
||
&+\underbrace{\langle \psi_i \vert \hat{T}_\mathrm{no}[\rho(\mathbf{r})]\vert \psi_i \rangle + \int \frac{\delta E_\mathrm{eff}[\rho(\mathbf{r})]}{\delta\rho(\mathbf{r})}\Big\vert _{\rho_0} \rho(\mathbf{r})d\mathbf{r}}_{E_\mathrm{band}} \nonumber \\
|
||
&+ \underbrace{\frac{1}{2} \iint \frac{\delta^2 E_\mathrm{eff}[\rho(\mathbf{r})]}{\delta\rho(\mathbf{r})\delta\rho(\mathbf{r'})} \Big\vert _{\rho_{_0}, ~ \rho_{_0}'} \delta\rho(\mathbf{r}) \delta\rho(\mathbf{r'})d\mathbf{r}d\mathbf{r'}}_{E_\mathrm{2nd}}
|
||
\label{Edftb-long}
|
||
\end{align}
|
||
|
||
From the first line of equation~\ref{Edftb-long}, we can define that the reference Hamiltonian $\hat{H}^0$ only depends on the reference
|
||
electron density $\rho_{_0}$:
|
||
\begin{eqnarray}
|
||
\label{H0}
|
||
\hat{H}^0=-\frac{1}{2}\nabla^2 + \underbrace{V_\mathrm{ext}[\rho(\mathbf{r})]+ \int \frac{\rho_{_0}'(\mathbf{r'})}{|\mathbf{r}-\mathbf{r'}|}d\mathbf{r'} +V_\mathrm{xc}[\rho_{_0}(\mathbf{r})] }_{V_\mathrm{eff}([\rho_{_0}(\mathbf{r})])}
|
||
\end{eqnarray}
|
||
where we combined the last three terms in the operator as $V_\mathrm{eff}([\rho_{_0}(\mathbf{r})])$.
|
||
|
||
The right-hand side terms in the first line of equation~\ref{Edftb-long} vary linearly with $\rho_0(\mathbf{r})$, and correspond to a
|
||
repulsive contribution $E_\mathrm{rep}$. The sum of the terms in the second line is the so-called band energy $E_\mathrm{band}$. The third
|
||
line is the second-order energy $E_\mathrm{2nd}$. Equation~\ref{Edftb-long} can be rewritten as follows:
|
||
\begin{align}
|
||
E_\mathrm{DFTB}[\rho_{_0}(\mathbf{r})+\delta \rho(\mathbf{r})]=& E_\mathrm{rep}[\rho_{_0}(\mathbf{r})] + \overbrace{\sum_{i}^\mathrm{occ} n_i \left\langle \psi_i \left\vert \hat{H}^0 \right\vert \psi_i \right\rangle}^{E_\mathrm{band}} ~ + E_\mathrm{2nd}\left[\rho_{_0}(\mathbf{r}), (\delta\rho(\mathbf{r}))^2\right]
|
||
\label{Eterms}
|
||
\end{align}
|
||
|
||
\textbf{Band energy term.}
|
||
In DFTB, one relies on the use of LCAO for the description of the Kohn-Sham molecular orbitals $\psi_i$(\textbf{r}). Here the atomic
|
||
orbitals are limited to the valence orbitals of atoms:
|
||
\begin{align}
|
||
\psi_i(\mathbf{r}) = \sum_\nu C_{i\nu} \varphi_\nu(\mathbf{r}-\mathbf{R_\alpha})
|
||
\label{orbitals}
|
||
\end{align}
|
||
$E_\mathrm{band}$ is the sum over the energies of all occupied orbitals obtained by diagonalization of the parameterized Hamiltonian matrix.
|
||
For atoms $\alpha$ and $\beta$, the corresponding atomic orbitals are $\varphi_\mu$ and $\varphi_\nu$.
|
||
Then $E_\mathrm{band}$ can be rewritten with equations~\ref{Veff} and \ref{orbitals}:
|
||
|
||
\begin{align}
|
||
\sum_{i}^\mathrm{occ} n_i \left\langle \psi_i \left\vert \hat{H}^0 \right\vert \psi_i \right\rangle=\sum_{i}^\mathrm{occ} n_i \sum_{\mu}^\mathrm{occ} \sum_{\nu}^\mathrm{occ} C_{i\mu}C_{i\nu} \underbrace{\left\langle \varphi_\mu \left\vert \hat{T}_e[\rho(\mathbf{r})] + V_\mathrm{eff}[\rho_{_0}(\mathbf{r})] \right\vert \varphi_\nu \right \rangle}_{H^0_{\mu\nu}}
|
||
\label{Eband}
|
||
\end{align}
|
||
$\forall\mu\in {\alpha}, \nu\in{\beta}$. The Hamiltonian matrix element $H^0_{\mu\nu}$ is defined as:
|
||
|
||
\begin{align}
|
||
H^0_{\mu\nu}&= \left\langle \varphi_\mu \left\vert \hat{H}^0 \right\vert \varphi_\nu \right\rangle, \nonumber \\
|
||
\label{matrix}
|
||
\end{align}
|
||
|
||
|
||
The effective potential $V_\mathrm{eff}[\rho_{_0}(\mathbf{r})]$ is defined as the sum of potentials $V_\alpha(\textbf{r})$ centered on the atoms:
|
||
\begin{align}
|
||
V_\mathrm{eff}[\rho_{_0}(\mathbf{r})]=\sum_{\alpha}V
|
||
_\alpha(\mathbf{r}_\alpha)
|
||
\label{Vdftb}
|
||
\end{align}
|
||
where $\mathbf{r}_\alpha=\mathbf{r}-\mathbf{R}_\alpha$.
|
||
The Hamiltonian matrix elements can be written as follows:
|
||
\begin{align}
|
||
H^0_{\mu\nu}=\left\langle\varphi_\mu \left\vert -\frac{1}{2}\nabla_\nu^2 + V_\alpha+(1-\delta_{\alpha\beta})V_\beta \right\vert \varphi_\nu \right\rangle, \mathrm {with} ~ \mu \in {\alpha}, \nu \in {\beta}
|
||
\label{Hamil}
|
||
\end{align}
|
||
where $\delta_{\alpha\beta}$ is the Kronecker's delta.
|
||
|
||
%In practice, $H^0_{\mu\nu}$ are calculated as follows:
|
||
For diagonal elements, the energy level in the free atom is chosen, which ensures the correct dissociation limits.
|
||
The interatomic blocks are computed as given in equation~\ref{Hamil}, depending on the choice of potential generation.
|
||
Because of the orthogonality of the basis functions, the off-diagonal elements of the intraatomic blocks are exactly zero.
|
||
To summarize, within the electronic density superposition approach, the $H^0_{\mu\nu}$ elements can be unfolded as:
|
||
\begin{align}
|
||
H^0_{\mu\nu}=
|
||
\begin{cases}
|
||
\varepsilon_\mu^\mathrm{free~atom}, ~ & \mu=\nu \\
|
||
\left \langle \varphi_\mu \left\vert -\frac{1}{2}\nabla_\nu^2 + V_\alpha+(1-\delta_{\alpha\beta})V_\beta \right\vert \varphi_\nu \right\rangle, ~ & \mu \in {\alpha}, \nu \in {\beta}, \alpha \neq \beta \\
|
||
0, ~ & \mathrm{otherwise}
|
||
\end{cases}
|
||
\label{Hmatrix}
|
||
\end{align}
|
||
|
||
It should be noted that the $H^0_{\mu\nu}$ elements only depend on atoms $\alpha$ and $\beta$. Therefore only the two-center matrix elements and
|
||
the two-center elements of the overlap matrix can be explicitly calculated, in other words, interactions at three or more centers are neglected as stated above.
|
||
|
||
\textbf{Repulsive energy term.} $E_\mathrm{rep}$ is a repulsive contribution obtained from the sum of atomic-pair terms, which only depend
|
||
on the reference electronic density $\rho_{_0}(\textbf{r})$.
|
||
|
||
%In order to obtain a good evaluation of $\rho_{_0}(\textbf{r})$, it is approximated as a overlap of atom-like densities centered on the nuclei $\alpha$:
|
||
%\begin{align}
|
||
%\rho_{_0}(\mathbf{r})=\sum_{\alpha}^{N}
|
||
%\rho_{_0}^\alpha(\mathbf{r}_\alpha), \mathrm{with} ~~
|
||
%\mathbf{r}_\alpha=\mathbf{r}-\mathbf{R}_\alpha
|
||
%\label{p0}
|
||
%\end{align}
|
||
%with this approximation, it ensures $E_\mathrm{rep}$ doesn't rely on the electronic density fluctuations.
|
||
%In addition, owing to the neutrality of $\rho_{_0}(\textbf{r})$ and Coulomb contributions can be negligible for long distances. So $E_\mathrm{rep}$ is expanded as follows:
|
||
%\begin{align}
|
||
%E_\mathrm{rep}=&\sum_{\alpha}^{N}E_\mathrm{rep}[
|
||
%\rho_{_0}^\alpha(\mathbf{r}_\alpha)]+\frac{1}{2}
|
||
%\sum_{\beta}^{N}\sum_{\alpha \neq \beta}^{N} \left(E_\mathrm{rep}[\rho_{_0}^\alpha(\mathbf{r}_\alpha)+\rho_{_0}^\beta(\mathbf{r}_\beta)]-E_\mathrm{rep}[\rho_{_0}^\alpha(\mathbf{r}_\alpha]-E_\mathrm{rep}[\rho_{_0}^\beta(\mathbf{r}_\beta] \right) \nonumber \\
|
||
%& +E_\mathrm{long}
|
||
%\label{E_{rep}}
|
||
%\end{align}
|
||
%where $E_\mathrm{long}$ refers to the interactions energy of more than two centers.
|
||
|
||
Due to the screening of terms of more than two centers, one can assume that the two-center contributions are short ranged.
|
||
But $E_\mathrm{rep}$ doesn't decay to zero for long interatomic distances. Instead, it decays to a constant given by the atomic
|
||
contributions:
|
||
\begin{align}
|
||
\displaystyle \lim_{R_{\alpha\beta} \to \infty}E_\mathrm{rep}[\rho_{_0}(\mathbf{r})]=\sum_{\alpha}^{N}
|
||
E_\mathrm{rep}[\rho_{_0}^\alpha(\mathbf{r}_\alpha)]
|
||
\label{Erep-infinity}
|
||
\end{align}
|
||
|
||
The right side of equation~\ref{Erep-infinity} is assumed to make $E_\mathrm{rep}$ only rely on the two-center contributions:
|
||
\begin{align}
|
||
E_\mathrm{rep}[\rho_{_0}(\mathbf{r})]\approx\frac{1}{2}\sum_{\alpha}^{N}\sum_{\beta}^{N} V(\mathbf{R}_\alpha-\mathbf{R}_\beta)
|
||
\label{Erep-right}
|
||
\end{align}
|
||
|
||
In practice, it it possible to calculate $E_\mathrm{rep}$ with known values of $\rho_{_0}(\textbf{r})$, but it's more convenient to adjust
|
||
the expression of $E_\mathrm{rep}$ to $ab ~ initio$ calculations. Therefore, $E_\mathrm{rep}$ is determined by comparing the
|
||
difference between the DFT energy $E_\mathrm{DFT}$ and $E_\mathrm{band}$+$E_\mathrm{2nd}$ as a function of the interatomic distance $R_{\alpha\beta}$:
|
||
\begin{align}
|
||
E_\mathrm{rep}[\rho_{_0}(\mathbf{r})] \equiv E_\mathrm{rep}(R_{\alpha\beta})=E_\mathrm{DFT}(R_{\alpha\beta})-E_\mathrm{band}(R_{\alpha\beta})-E_\mathrm{2nd}(R_{\alpha\beta})
|
||
\label{Erep-final}
|
||
\end{align}
|
||
|
||
\textbf{Second-order term.}
|
||
In SCC-DFTB, the electronic density is corrected by including the second-order contribution $E_\mathrm{2nd}$ in equation~\ref{Eterms},
|
||
which is ignored in first-order DFTB.
|
||
|
||
To include the density fluctuations in a simple but efficient way according to tight-binding method, $\delta \rho$ can be
|
||
written as a superposition of atom-like contributions $\delta \rho_\alpha$, which has a fast decrease with the increase of
|
||
the distance from the corresponding atomic center:
|
||
|
||
\begin{align}
|
||
\delta\rho(\mathbf{r})=\sum_{\alpha}^{N}
|
||
\delta\rho_\alpha(\mathbf{r})
|
||
\label{relative_p}
|
||
\end{align}
|
||
where $\delta \rho_\alpha$ can be simplified with the monopole approximation as follows:
|
||
|
||
\begin{align}
|
||
\delta\rho_\alpha(\mathbf{r})=\Delta q_\alpha F_{0}(\mathbf{r}-\mathbf{R}_\alpha)
|
||
\label{relative_pa}
|
||
\end{align}
|
||
where the atomic charge fluctuation $\Delta q_\alpha$ (difference between the Mulliken population $q_\alpha$ \cite{Mulliken1955} of atomic $\alpha$ and the number of valence electrons of the atom at infinity) is estimated by the Mulliken expression. $F_{0}$ represents the normalized radial dependence of the electronic density fluctuation in atom $\alpha$. %approximated to spherical by the angular function $Y_{00}$.
|
||
This means the effects of charge transfer are included, however, the changes in the shape of the electronic density are ignored.
|
||
$E_\mathrm{2nd}$ can be rewritten with equations~\ref{Edftb-long} and \ref{Eterms} as follows:
|
||
|
||
\begin{align}
|
||
E_\mathrm{2nd}\approx
|
||
&\frac{1}{2}\sum_{\alpha}^{N}\sum_{\beta}^{N} \Delta q_{\alpha} \Delta q_{\beta} \overbrace{\iint \left(\frac{1}{|\mathbf{r}-\mathbf {r'}|} + \frac{\delta^2 E_{xc}[\rho_{_0}(\mathbf{r})]}{\delta\rho(\mathbf{r})\delta\rho(\mathbf{r'})} \Big\vert _{\rho_{_0}, ~ \rho_{_0}'} \right) F(\alpha, \beta) d\mathbf{r}d\mathbf{r'}}^{\gamma_{\alpha\beta}} \nonumber \\
|
||
=& \frac{1}{2}\sum_{\alpha}^{N}\sum_{\beta}^{N} \Delta q_{\alpha} \Delta q_{\beta} \gamma_{\alpha\beta} \\
|
||
F(\alpha, \beta)= & F_{0}(\mathbf{r}-\mathbf{R}_\alpha)\times F_{0}(\mathbf{r'}-\mathbf{R}_\beta)
|
||
\label{E2nd}
|
||
\end{align}
|
||
where the two-electron integrals $\gamma_{\alpha\beta}$ is introduced for convenience.
|
||
|
||
To calculate equation~\ref{E2nd}, $\gamma_{\alpha\beta}$ must be analyzed. In the limit case, the interatomic distance is very large,
|
||
$|{\textbf{R}_\alpha - \textbf{R}_\beta}|$ = $|\textbf{r} - \textbf{r}'| \rightarrow \infty$ with GGA-DFT, the exchange-correlation term tends to zero.
|
||
$\gamma_{\alpha\beta}$ that describes the interaction of two normalized spherical electronic densities reduces to 1/$|\textbf{R}_\alpha-\textbf{R}_\beta|$, so $E_\mathrm{2nd}$ can be expressed as follows:
|
||
|
||
\begin{align}
|
||
E_\mathrm{2nd}=\frac{1}{2}\sum_{\alpha}^{N}\sum_{\beta}^{N} \frac{\Delta q_{\alpha} \Delta q_{\beta}} {|\mathbf{R}_\alpha-\mathbf{R}_\beta|}
|
||
\label{E2ndsimple}
|
||
\end{align}
|
||
|
||
It is worth noting that the electronic density $\rho(\mathbf{r})$ influences explicitly the calculation of the electrostatic energy in DFT.
|
||
In the context of DFTB, point charges are used and the electronic density around the atom is condensed at a point. In practice, Mulliken's
|
||
definition of charge is often used,\cite{Mulliken1955} which is defined as:
|
||
\begin{align}
|
||
q_\alpha^\mathrm{Mull}=\frac{1}{2}\sum_{\alpha}^\mathrm{occ} n_i
|
||
\sum_{\mu \in \alpha} \sum_{\nu} \left(C_{i\mu}^*C_{i\nu}S_{\mu\nu} + C_{i\nu}^*C_{i\mu}S_{\nu\mu}
|
||
\right)
|
||
\label{Mulliken}
|
||
\end{align}
|
||
where $\mu$ denotes the orbitals belonging to the atom $\alpha$. This definition does not allow the bond between two different atoms to polarize.
|
||
|
||
\textbf{Total energy.}
|
||
The total energy in SCC-DFTB can be written from the previous different contributions as follows:
|
||
\begin{align}
|
||
E_\mathrm{SCC}=
|
||
& \sum_{i}^\mathrm{occ} n_i \sum_{\mu}^\mathrm{occ} \sum_{\nu}^\mathrm{occ} C_{i\mu}C_{i\nu} \left\langle \psi_i \left\vert \hat{H}^0 \right\vert \psi_i \right\rangle + \frac{1}{2}\sum_{\alpha}^{N}\sum_{\beta}^{N}
|
||
V(\mathbf{R}_\alpha-\mathbf{R}_\beta) \nonumber \\
|
||
& + \frac{1}{2}\sum_{\alpha}^{N}\sum_{\beta}^{N} \Delta q_{\alpha} \Delta q_{\beta} \gamma_{\alpha\beta}
|
||
\label{Escc}
|
||
\end{align}
|
||
%where $\gamma_{\alpha\beta}$=$\gamma_{\alpha\beta}(V_\alpha, V_\beta, \mid{\textbf{R}_\alpha - \textbf{R}_\beta}\mid)$.
|
||
%<<<<<<< HEAD
|
||
%Here the contribution from $H^0$ is exactly the same with the one in standard DFTB scheme and the wavefunctions $\psi_i$ are expanded in a LCAO model (eq \ref{orbitals}).
|
||
|
||
%\textbf{Secular equation}
|
||
%=======
|
||
Here the contribution from $H^0$ is exactly the same with the one in standard DFTB scheme and the molecular orbitals $\psi_i$ are expanded in a LCAO model (equation~\ref{orbitals}).
|
||
%>>>>>>> 92023a10c3aa8b7dc4ace43987c1d571fb99a738
|
||
|
||
\textbf{Secular equations.}
|
||
From this LCAO model, we can get the secular equations:
|
||
\begin{align}
|
||
\sum_\nu C_{i\nu} \left({H}^0_{\mu\nu} - \varepsilon_i S_{\mu\nu} \right) =0, ~ \forall \mu, ~ \nu
|
||
\label{secular-eq}
|
||
\end{align}
|
||
where $H^0_{\mu\nu}$ are the Hamiltonian matrix elements and $S_{\mu\nu}$ are the overlap matrix elements.
|
||
|
||
The secular equations can be rewritten with modified Hamiltonian matrix elements, $H_{\mu\nu}$, defined by:
|
||
\begin{align}
|
||
H_{\mu\nu}=H_{\mu\nu}^0 + \frac{1}{2} S_{\mu\nu} \sum_{\zeta}(\gamma_{\alpha\zeta} + \gamma_{\beta\zeta})\Delta q_\zeta
|
||
\label{Huv}
|
||
\end{align}
|
||
|
||
The $H_{\mu\nu}^0$ and $S_{\mu\nu}$ matrix elements are the same than the ones defined in the standard DFTB method (equation~\ref{matrix}).
|
||
The $H_{\mu\nu}$ elements rely on the atomic charges explicitly, and the atomic charges depend on the molecular orbitals (see equation~\ref{Mulliken}).
|
||
Then the resolution can be achieved in a self-consistent way.
|
||
First, from an initial set of charges
|
||
%$\lbrace q_\alpha \rbrace$,
|
||
the $H_{\mu\nu}$ elements which depend on these charges can be computed.
|
||
The KS equations (equations~\ref{secular-eq}) are then solved which gives the energy of the KS orbitals and the corresponding eigenvectors.
|
||
The corresponding coefficients allow to compute a new set of charges which will be used in the calculation of new $H_{\mu\nu}$ elements.
|
||
This procedure is repeated until the atomic charges are converged.
|
||
%The self-consistent charge correction allows to treat the charge transfer effects explicitly in which the transfer-ability of $E_{rep}$ is fairly better compared to the non-self-consistent scheme.
|
||
|
||
DFTB is derived from DFT, it therefore inherits the specific problems of DFT. For instance, the traditional DFT functionals can not describe properly
|
||
dispersion interaction and charge resonance phenomena in charged aggregates. DFTB also display some specific problems because of its own
|
||
approximations such as the use of Mulliken atomic charges, the absence of atomic polarization, the absence of coupling between atomic orbitals
|
||
located on the same atom. This differs from DFT that explicitly considers atomic polarization.
|
||
|
||
\textbf{Atomic charges.}
|
||
As presented above, DFTB was initially developed with Mulliken charges; However, other definitions of atomic charges are possible such as
|
||
Natural Bond Order (NBO)\cite{Bader1985, Bader1990, Glendening1998, Glendening2012} and Electrostatic Potential Fitting (EPF) charges.\cite{Singh1984, Besler1990}
|
||
EPF has a fairly good representation of the electrostatic term of a molecule dominated by the Van der Waals interactions. CM3 (Class IV / Charge Model 3)
|
||
charges were proposed by Li $et al.$ in 1998 and they have been considered in DFTB. They give good results for the description of the electric dipole
|
||
and the electrostatic potential, partial atomic charges in molecules, and Coulombic intermolecular potential of polycyclic aromatic hydrocarbon clusters.\cite{Li1998, Kalinowski2004, Kelly2005, Rapacioli2009corr}
|
||
The CM3 charges are defined as:
|
||
\begin{align}
|
||
q_{\alpha}^\mathrm{CM3} = q_{\alpha}^\mathrm{Mull} + \sum_{\beta \neq \alpha} (C_{t_\alpha t_\beta} K_{\alpha\beta} + D_{t_\alpha t_\beta} K_{\alpha\beta}^2 )
|
||
\label{CM3}
|
||
\end{align}
|
||
where $K_{\alpha\beta}$ is the Mayer bond order,\cite{Mayer1983, Cances1997, Lu2012} which depends on the density matrix of the orbitals belonging to
|
||
each of the two atoms $\alpha$ and $\beta$. $C_{t_\alpha t\beta}$ and $D_{t_\alpha t\beta}$ are empirical parameters which are related
|
||
to the nature of atoms $\alpha$ and $\beta$.
|
||
|
||
In this thesis, in practice for the calculation of electronic energy, the definition of CM3 charges is simplified as:
|
||
\begin{align}
|
||
q_{\alpha}^\mathrm{CM3} = q_{\alpha}^\mathrm{Mull} + \sum_{\beta \neq \alpha} C_{t_\alpha t_\beta} K_{\alpha\beta}
|
||
\label{CM3-classical}
|
||
\end{align}
|
||
|
||
|
||
%<<<<<<< HEAD
|
||
%\textbf{Dispersion energy}
|
||
%
|
||
%In order to correctly describe the energies of molecular systems, it is necessary to take into account the Van der Waals forces. London dispersal force is a type of force acting between atoms and molecules
|
||
%that are normally electrically symmetric (symmetrically distributed electrons with respect to the nucleus), which is a large part of the Van der Waals force. London dispersion force arises from the interactions between fluctuating dipoles.
|
||
%One of the major problems within DFTB is that it does not take the dispersion interactions into account. The semi-empirical energy corrections have been applied on DFT computations, which usually gives good results. \cite{Lewis1995, Gianturco1999, Elstner2001, Wu2002, Grimme2004, Zimmerli2004, Neumann2005, Goursot2007}
|
||
%Moreover, a semi-empirical correction of the dispersion to the DFTB energy leaves the freedom to use the already existing DFTB parameters.
|
||
%For the study of water clusters and PAH dimers in this thesis, the dispersion energy plays a fairly important role.
|
||
%The correction of the dispersion energy is of the form as follows:
|
||
%=======
|
||
\textbf{Dispersion energy.}
|
||
In order to correctly describe the energies of molecular systems, it is necessary to take into account the Van der Waals interactions.
|
||
London dispersion interaction acts between atoms and molecules and represents a large part of the Van der Waals interactions.
|
||
London dispersion interaction arises from the interactions between fluctuating dipoles. One of the major drawbacks within DFTB
|
||
is that it does not take dispersion interactions into account. This is also true for DFT, when using LDA or GGA functionals.
|
||
To overcome this limitation, semi-empirical energy corrections can be applied in DFT and DFTB calculations, which usually gives
|
||
good results.\cite{Lewis1995, Gianturco1999, Elstner2001, Wu2002, Grimme2004, Zimmerli2004, Neumann2005, Goursot2007}
|
||
Moreover, a semi-empirical correction to the dispersion in the DFTB energy leaves the freedom to use the already existing DFTB
|
||
parameters withjout any need of re-parametrization. For the studies presented in this thesis, the dispersion energy plays a fairly
|
||
important role. The correction we apply for the dispersion energy is of the following form:
|
||
%>>>>>>> 92023a10c3aa8b7dc4ace43987c1d571fb99a738
|
||
%\cite{Kaplan2006, Koskinen2009}
|
||
\begin{align}
|
||
E_\mathrm{disp}(\mathbf R_{\alpha\beta})=-\sum_{\alpha} \sum_{\beta \neq \alpha} f(\mathbf R_{\alpha\beta})
|
||
\frac{C_{\alpha\beta}^6}{\mathbf R_{\alpha\beta}^6}
|
||
\label{dispersionE}
|
||
\end{align}
|
||
where $f(\textbf R_{\alpha\beta})$ is a cutoff function, which allows to avoid the divergence of this term at a short distance. $C_{\alpha\beta}^6$ is an empirical coefficient calculated for each pair of atoms.
|
||
|
||
|
||
\subsection{Force Field Methods}
|
||
|
||
Force field is a computational method utilized to estimate the forces between particles, in other words it is the functional form and parameter sets applied to calculate the potential energy of a system.
|
||
FF is interatomic potential and use the same concept with force field in classical physics, a vector field which describes a non-contact force acting on a particle at different positions in space.
|
||
The acting force on each particle is derived as a gradient of the potential energy with respect to the particle positions.\cite{Frenkel2001}
|
||
In such case, the interactions in a system are determined from parameterized potentials in which the electronic structure can not be described explicitly because each particle is treated as a material point. The particles interact with each other through the FF and the integration algorithm is applied to the particles. In most cases, this leads to a big decrease of the precision level in the description of the system but it can reduce the calculation cost drastically, which allows to model systems containing several thousands of particles.
|
||
|
||
Different potentials have been proposed,\cite{Jones1924, Lennard1924, Morse1929, Born1932, Stillinger1985, Kaplan2006} which can be classified two main groups: pair potentials and multi-body potentials. For pair potentials, harmonic interaction is the most basic form:\cite{Bahar1997}
|
||
\begin{align}
|
||
V(R_{\alpha\beta}) = k(R_{\alpha\beta}-R_{eq})^2
|
||
\label{harmonic}
|
||
\end{align}
|
||
where $R_{\alpha \beta}$ is the distance between two interacting particles $\alpha$ and $\beta$. $k$ is the harmonic force constant. $R_{eq}$ is the equilibrium distance where the force of repulsion equals to the one of attraction. This potential is sufficient for systems only deviating very slightly from the bond distance at equilibrium and interactions reasonably limited to adjacent pairs of particles. However, for systems with large deviations, other potential forms must be used, for instance, the Morse potential which describes the potential energy of a diatomic molecule:\cite{Morse1929, Girifalco1959}
|
||
\begin{align}
|
||
V(R_{\alpha\beta}) = D_{eq}\left(1-e^{-a(R_{\alpha\beta}-R_{eq})} \right)^2
|
||
\label{Morse}
|
||
\end{align}
|
||
where $D_{eq}$ is the depth of the Morse potential well. Parameter $a$ determines the width of the potential, the smaller $a$ the larger the well. The force constant of the bond can be found via the Taylor expansion of $V(R_{\alpha\beta})$ around $R_{\alpha\beta}=R_{eq}$ to the second derivative of the potential energy function, from which it can obtain
|
||
$a=(k_{eq}/{2D_{eq}})^{\frac{1}{2}}$ in which $k_{eq}$ is the force constant of the minimum well.
|
||
|
||
Lennard-Jones potential as known as LJ potential or 12-6 potential is a pair potential, which is proposed by John Lennard-Jones in 1924.\cite{Jones1924, Lennard1924} It models soft repulsive and attractive interactions, therefore, the LJ potential describes electronically neutral atoms or molecules.
|
||
Because of its simple mathematical form, it is one of the most widely used intermolecular potentials especially to describe the interaction within noble gas molecules. The total energy can be written as the sum of the interaction energy of all atomic pairs, which is defined as follows:\cite{Lennard1931}
|
||
\begin{align}
|
||
V_\mathrm{LJ}(R_{\alpha\beta}) = 4\varepsilon_0 \left[ \left(\frac{\sigma}{R_{\alpha\beta}}\right)^{12}- \left(\frac{\sigma}{R_{\alpha\beta}}\right)^6
|
||
\right] = \varepsilon_0 \left[ \left(\frac{R_{eq}}{R_{\alpha\beta}}\right)^{2n}- 2\left(\frac{R_{eq}}{R_{\alpha\beta}}\right)^n \right]
|
||
\label{LennarsJones}
|
||
\end{align}
|
||
$\varepsilon_0$ denotes the depth of the potential well that usually refers to the dispersion energy.
|
||
$\sigma$ is the interparticle distance at which the potential energy is zero.
|
||
$n = 6$ and $\varepsilon_0$ denotes the bonding energy, the energy required to separate the atoms. And the LJ potential has its minimum ($-\varepsilon_0$) at a distance of $R_{\alpha \beta}=R_{eq} = 2^{\frac{1}{6}}\sigma$.
|
||
|
||
It becomes more complicated for molecular systems in which the modes of intermolecular and intramolecular interactions are very different. So it requires to develop force fields including several kinds of potentials.
|
||
The expression of the potential energy for a molecular system which is the most frequently used for simple organic molecules and biological macromolecules is written as follows:\cite{Monticelli2013}
|
||
\begin{align}
|
||
V_\mathrm{total}(R) =& \overbrace{V_\mathrm{bond} + V_\mathrm{angle} + V_\mathrm{dihedral}}^{V_\mathrm{bonded}} + \overbrace{V_\mathrm{VW} + V_\mathrm{Coulomb}}^{V_\mathrm{nonbonded}} \nonumber \\
|
||
=& \sum_\mathrm{bond} \frac{k_{\alpha\beta}}{2}(R_{\alpha\beta}-R_{eq})^2 + \sum_\mathrm{angle} \frac{k_\theta}{2}(\theta-\theta_{eq})^2 + \sum_\mathrm{dihedral} \frac{k_\phi}{2}(1+ \mathrm{cos}(n \phi-\phi_{eq}))^2 \nonumber \\
|
||
& + \sum_\mathrm{VW} 4\varepsilon_0 \left[ \left(\frac{\sigma}{R_{\alpha\beta}}\right)^{12}- \left(\frac{\sigma}{R_{\alpha\beta}}\right)^6 \right] +
|
||
\sum_\mathrm{Coulomb} \frac{1}{4\pi \varepsilon_0}\frac{q_{\alpha}q_{\beta}}{R_{\alpha\beta}}
|
||
\label{potential}
|
||
\end{align}
|
||
Molecular interactions can determine the macroscopic properties of matter. Van der Waals interaction is an important force between atoms and molecules but it is extremely short ranged. Van der Waals interaction energy is also termed London dispersion energy.
|
||
This repulsive distance dependence is usually modeled as a term that scales with $1/R_{\alpha\beta}^{12}$ although there is no absolute physical reason for it.
|
||
The van der Waals interaction energy is usually approximated by a Lennard-Jones potential.
|
||
|
||
%The partial derivative of $V_\mathrm{LJ}(R_{\alpha\beta})$ with respect to the distance $R_{\alpha \beta}$ is:
|
||
%\begin{align}
|
||
%\frac{\partial V_\mathrm{LJ}(R_{\alpha\beta})}{\partial R_{\alpha \beta}} = -\frac{24 \varepsilon_0}{R_{\alpha \beta}} \left[ 2 \left(\frac{\sigma}{R_{\alpha\beta}}\right)^{12} - \left(\frac{\sigma}{R_{\alpha\beta}}\right)^6 \right]
|
||
%\label{partialLJ}
|
||
%\end{align}
|
||
%In this case, the force acting on the atom $\alpha$ takes the following form:
|
||
%\begin{align}
|
||
%\mathbf{F}_\alpha^\mathrm{VW} = \frac{24 \varepsilon_0}{R_{\alpha \beta}} \left[2 \left(\frac{\sigma}{R_{\alpha\beta}}\right)^{12} - \left(\frac{\sigma}{R_{\alpha\beta}}\right)^6 \right]\mathbf{u}=-\mathbf{F}_\beta^\mathrm{VW}
|
||
%\label{VanderWaF}
|
||
%\end{align}
|
||
%where $\textbf{u}$ is the normalization of $\overrightarrow{\alpha\beta}$. From the action and reaction principle, we can deduce the force on atom $\beta$ is the opposite direction of that on atom $\beta$.
|
||
|
||
Now we know how to calculate the potential energy of nuclei (electronic energy).
|
||
Many thermal dynamical, chemical, and physical properties require the association of PES.
|
||
Different methods will be shown for the exploration of the PES in the next section.
|
||
|
||
%<<<<<<< HEAD
|
||
%\section{Exploration of the PES}
|
||
|
||
%The PES is a function giving the energy of a system, according to one or more nuclear coordinates. If there is only one coordinate, the PES is called the energy profile or potential energy curve. Figure \ref{PES} is a model of PES.
|
||
%The analogy of a hilly landscape with peaks, valleys, mountainpasses helps to understand the PES. Stable molecular structures correspond to the positions of the minima in the valleys on a PES. The shape of the valley around a minimum determines the vibrational spectrum.
|
||
%The points on a PES can be classified according to the first and all second derivatives of the energy with respect to position, which are the gradient and the curvature, respectively.
|
||
%When points with a zero gradient (stationary points) and the second derivatives are positive, this corresponds to a local minimum (physically stable structure). Among these local minima, the lowest energy minima is called the global minimum. When at least one of the second derivatives is negative, it corresponds to a transition state (saddle point), which is showed in Figure \ref{PES}.
|
||
|
||
%\figuremacro{PES}{Some key point representations on a PES.}
|
||
|
||
%The molecular structure, properties, chemical reactivity dynamics, and spectra of molecules can be readily understood in terms of PES. The concept of PES can be applied in many fields, for instance, the physics and chemistry especially in the theoretical subbranches of these subjects. Only very simple PES can be obtained from experiment while the computational chemistry field has developed different kinds of methods to explore the PES.
|
||
%To survey the PES, the choice of the exploration methods can be guided by the form of the PES (the statistical set that one wishes to study) and the temporal aspect.
|
||
%%Marcella $et al.$ proposed that Car-Parrinello MD can do an efficient exploration for reactive PES.\cite{Iannuzzi2003}
|
||
%%The PES was also explored by a machine learning method to characterize the atomic transport. \cite{Kanamori2018}
|
||
|
||
%Monte Carlo methods and MD methods are widely used efficient approaches for the exploration of PES of systems containing a large number of degrees of freedom such as molecular aggregates. Monte Carlo methods allow the sampling of a PES by performing random shifts in order to correctly reproduce the probability distribution of the configurations which are accessible from the phase space of a system as a function of state variables such as temperature, energy, or number of particles. One of the main advantages of the Monte Carlo methods is that it does not require the calculation of gradients, so it avoids a major drawback of MD algorithms. In MD simulations, it needs to numerically solving Newton’s equations of the interacting many-body system, and one can obtain static properties by taking averages along the resulting deterministic trajectory in phase space. Parallel-tempering molecular dynamics (PTMD) is an improvement of classical MD. The principles of Monte Carlo, classical MD and PTMD methods will be described.
|
||
|
||
%\subsection{Monte Carlo Simulations}
|
||
|
||
%The term Monte Carlo denotes a class of algorithmic methods that aim at a probabilistic description, relying on the use of random numbers. The name of Monte Carlo alludes to the games of chance taking place at the Monte Carlo casino. The Monte Carlo method was introduced in 1947 by Metropolis,\cite{Metropolis1987} and was first published in 1949 by Metropolis collaborated with Ulam.\cite{Metropolis1949 Monte Carlo methods have been widely applied in computational physics, computational statistics, biomedicine, machine learning, industrial engineering, economics and finance, and other fields.\cite{Kroese2014} Monte Carlo methods can generally be roughly divided into two categories. One is that the problem to be solved has inherent randomness, and the computing power of the computer can directly simulate this random process. For example, in the research of nuclear physics, analysis of the transmission process of neutrons in the reactor.
|
||
%The other type is that the problem to be solved can be transformed into a randomly distributed characteristic numbers, such as the probability of a random event. Through the random sampling method, the probability of random events is estimated by the frequency of occurrence, and this is used as the solution to the problem. This method is usually used to solve complicated multi-dimensional integration problems.
|
||
%
|
||
%Monte Carlo methods are mainly used in the problems of optimization, numerical integration, and generating draws from a probability distribution.
|
||
%In physics related problems Monte Carlo methods are applied in the simulated systems with many coupled degrees of freedom, such as fluids, strongly coupled solids, disordered materials, and cellular structures, which are useful when it is hard or impossible to use other approaches. In statistical physics not related to thermodynamics, Monte Carlo molecular simulation is an alternative to MD simulation, and Monte Carlo methods can be used to compute statistical field theories of simple particle and polymer systems.\cite{Rosenbluth1955, Binder1993, Baeurle2009}
|
||
%Monte Carlo simulations are the generation of random objects or processes.
|
||
%When solving practical problems using Monte Carlo method, it includes two main parts. First, random variables with various probability distributions need to be generated to simulate a certain process with Monte Carlo method. Second, statistical methods need to be used to estimate the numerical characteristics of the model to obtain a numerical solutions to the actual problem.
|
||
%Monte Carlo methods make it possible for the sampling of a PES by performing random shifts in order to correctly reproduce the probability distribution of the configurations.
|
||
%One example for the calculation procedures of molecular simulation are as follows:
|
||
|
||
%1. Initially, a random molecular configuration is generated using a random number generator.
|
||
%=======
|
||
\section{Exploration of PES}
|
||
|
||
The PES is a function giving the energy of a system according to one or more nuclear coordinates.
|
||
If there is only one coordinate, the PES is called the energy profile or potential energy curve. Figure~\ref{PES}
|
||
is a model of PES in two dimension. The analogy of a hilly landscape with peaks, valleys, mountain passes helps
|
||
to understand the PES. Stable molecular structures correspond to the minima in the valleys on a PES. The shape
|
||
of the valley around a minimum determines the vibrational spectrum.
|
||
The key points on a PES can be classified according to the first and all second derivatives of the energy with
|
||
respect nuclear coordinates, which correspond to the gradients and the curvatures, respectively.
|
||
When points have a zero gradient (stationary points) and their second derivatives are positive, this corresponds
|
||
to a \textbf{local minimum} (physically stable structure). Among these local minima, the lowest energy minima is called
|
||
the \textbf{global minimum}. When at least one of the second derivatives is negative, the point is a transition state (saddle point).
|
||
|
||
\figuremacro{PES}{Schematic representation of some key points on a model potential energy surface.}
|
||
|
||
Molecular structure, properties, chemical reactivity, dynamics, and vibrational spectra of molecules can be readily
|
||
understood in terms of PES. Only very simple PES can be obtained from experiment whereas computational chemistry
|
||
has developed different kinds of methods to efficiently explore PES. To survey the PES, the choice of an exploration
|
||
method can be guided by the shape of the PES (the statistical set that one wishes to study) and the temporal aspect.
|
||
|
||
\textbf{Monte Carlo} and classical \textbf{MD} simulations are widely recognized approaches for the exploration of PES of systems containing a large number of degrees of freedom such as molecular aggregates.
|
||
Monte Carlo methods allow the sampling of a PES by performing random shifts in order to correctly reproduce the
|
||
probability distribution of the configurations which are accessible from the phase space of a system as a function
|
||
of state variables such as temperature, energy, or number of particles.
|
||
One of the main advantages of Monte Carlo methods is that they do not require the calculation of gradients.
|
||
In MD simulations, one needs to numerically integrate Newton’s equations of motion of the interacting particles.
|
||
One can then obtain statistical properties by performing time averaging along the resulting deterministic trajectory
|
||
in phase space. \textbf{Parallel-tempering molecular dynamics} (PTMD) is an improvement of classical MD that allows
|
||
for a more efficient sampling of the PES. The Monte Carlo, classical MD and PTMD methods are briefly described in the following sections.
|
||
|
||
\subsection{Monte Carlo Simulations}
|
||
|
||
The term Monte Carlo denotes a class of algorithmic methods that aims at a probabilistic description, relying on the use of
|
||
random numbers. The name Monte Carlo alludes to the games of chance taking place at the Monte Carlo casino. The
|
||
Monte Carlo method was introduced in 1947 by Metropolis,\cite{Metropolis1987} and was first published in 1949 by
|
||
Metropolis in collaboration with Ulam.\cite{Metropolis1949} Monte Carlo methods have been widely applied in computational
|
||
physics, computational statistics, biomedicine, machine learning, industrial engineering, economics and finance, and other fields.\cite{Kroese2014}
|
||
\textbf{Monte Carlo methods can} generally be roughly \textbf{divided into two categories}. The first category is applied to problems that have
|
||
inherent randomness, and the computing power of the computer can directly simulate this random process. For example,
|
||
in nuclear physics, analysis of the transmission process of neutrons in a reactor.
|
||
The second category applies to problems that can be transformed into randomly distributed characteristic numbers, such as the
|
||
probability of a random event. Through the random sampling method, the probability of random events is estimated by the
|
||
frequency of occurrence, and this is used as the solution to the problem. This method is usually used to solve complicated
|
||
multi-dimensional integration problems.
|
||
|
||
Monte Carlo methods are mainly used in the field of optimization, numerical integration, and generating draws from a probability
|
||
distribution. In physics related problems, Monte Carlo methods are applied to simulat systems diplaying many coupled degrees
|
||
of freedom, such as fluids, strongly coupled solids, disordered materials, and cellular structures. In statistical physics not
|
||
related to thermodynamics, Monte Carlo molecular simulation is an alternative to MD simulations. Monte Carlo methods
|
||
can be used to compute statistical field theories of simple particles and polymer systems.\cite{Rosenbluth1955, Binder1993, Baeurle2009}
|
||
When solving practical problems using a Monte Carlo methods, it includes two main parts. First, random variables with various
|
||
probability distributions need to be generated to simulate a certain process. Second, statistical methods need to be used to
|
||
estimate the numerical characteristics of the model to obtain a numerical solutions to the actual problem. Monte Carlo methods
|
||
make it possible for the sampling of a PES by performing random shifts in order to correctly reproduce the probability distribution
|
||
of the configurations. One example for the calculation procedure of a molecular simulation is as follows:
|
||
|
||
\begin{enumerate}
|
||
\item A random molecular configuration is generated using a random number generator.
|
||
|
||
\item Random changes are made to the particle coordinates of this molecular configuration, resulting in a new molecular configuration.
|
||
|
||
\item Calculate the energy of the new molecular configuration.
|
||
|
||
\item Compare the energy change between the new molecular configuration and the former one to determine whether to accept the configuration change or not.
|
||
%>>>>>>> 92023a10c3aa8b7dc4ace43987c1d571fb99a738
|
||
|
||
\begin{itemize}
|
||
\item[$\bullet$] If the energy of the new molecular configuration is lower than that of the original one, the new configuration change is accepted, and the
|
||
new configuration is used for the next iteration.
|
||
|
||
\item[$\bullet$] If not, the Boltzmann factor is calculated and a random number is generated. If this random number is greater than the calculated
|
||
Boltzmann factor, then the new configuration is discarded. If not, the new configuration is conserved and is used for the next iteration.
|
||
|
||
\end{itemize}
|
||
|
||
\item If a new iteration is required, the process is repeated from step 2.
|
||
|
||
\end{enumerate}
|
||
|
||
|
||
\textbf{General principles.}
|
||
Monte Carlo calculations that lead to quantitative results may be regarded as attempts to estimate the value of a multiple integral.
|
||
This is particularly true for the applications in equilibrium statistical thermodynamics, where one hopes to calculate the thermal
|
||
average $\langle A \rangle _T$ of an observable $A(\textbf X)$ as an integral over phase space $\Omega$, where $\textbf X$ is a
|
||
point in $\Omega$:
|
||
\begin{align}
|
||
\langle A \rangle _T= \frac{1}{Z} \int_\Omega d\mathbf X A(\mathbf X)e^{-H(\mathbf X)/k_BT}
|
||
\label{thermalAverage}
|
||
\end{align}
|
||
in which $Z$ is the partition function, and $k_B$ is the Boltzmann constant. $T$ refers to the temperature, and $H(\textbf X)$ denotes the Hamiltonian of the system.
|
||
To illustrate the general application of Monte Carlo techniques, we take here the standard example of the one-dimensional integral $I$ over integration space $\Omega$:
|
||
\begin{align}
|
||
I=\int_\Omega f(\mathbf x)d\mathbf x
|
||
\label{SingleIntegral}
|
||
\end{align}
|
||
|
||
Such integral can be evaluated more efficiently using conventional numerical means than using Monte Carlo methods. However, the extension to higher
|
||
dimensions (number of degrees of freedom greater than 3) is always difficult in practice with conventional numerical integration. It is then possible to
|
||
use stochastic approaches that one explores the configuration space randomly and the computation of the integral is estimated in the form of an average
|
||
value:
|
||
\begin{align}
|
||
I_N=\frac{V}{N} \sum_i^N f(\mathbf x_i)
|
||
\label{SingleAverage}
|
||
\end{align}
|
||
where $V$ is the volume of the integration space and $N$ is the number of points drawn randomly.
|
||
Within the limit of numbers ($N \rightarrow \infty$), the computation of the integral is exact as follows:
|
||
\begin{align}
|
||
I= \displaystyle \lim_{N \rightarrow \infty} I_N
|
||
\label{limI}
|
||
\end{align}
|
||
In addition, according to the central limit theorem, we can obtain a Gaussian distribution of the integral result and the statistical error is independent of the dimension of the problem to be solved.
|
||
|
||
In practice, a completely random exploration is inefficient if $f(\mathbf x)$ is located in a specific space region. It is then interesting to constrain the exploration of space by introducing a distribution $\rho(\mathbf x)$ to preferentially visit the regions where the function has the maximum influence in the integral. The integral can be rewritten as:
|
||
\begin{align}
|
||
I=\int_\Omega \frac{f(\mathbf x)}{\rho(\mathbf x)} \rho(\mathbf x)d\mathbf x
|
||
\label{IntegralPx}
|
||
\end{align}
|
||
then the integral is estimated with the mean value of $f/p$ of the points explored:
|
||
\begin{align}
|
||
I_N=\frac{V_p}{N} \sum_i^N \frac{f(\mathbf x_i)}{\rho(\mathbf x_i)}
|
||
\label{INp}
|
||
\end{align}
|
||
where $V_p$ is the volume of the integration space weighted by the distribution $\rho(\mathbf x)$. In order to explore the surface as efficiently as possible in calculating the mean value, the best choice for the points distribution is $\rho(\mathbf x) = f(\mathbf x)$.
|
||
|
||
In physical system, it usually needs to consider the weighting of different possible configurations to calculate the mean values.
|
||
By definition, the computation of an average value of an observable $A$ is the integral of $A$ over the whole phase space weighted by $\rho(\mathbf x)$ divided by the volume of this weighted space and is written as follows:
|
||
\begin{align}
|
||
\langle A \rangle= \frac{\int A(\mathbf x)\rho(\mathbf x) d\mathbf x}{\int \rho(\mathbf x)d \mathbf x}
|
||
\label{averaA}
|
||
\end{align}
|
||
Then the problem is to generate a distribution of configurations according to the law $\rho(\mathbf x)$.
|
||
The solution used by the Monte Carlo methods is to generate these points using a Markov chain, in other words, a sequential series of configurations where each configuration belongs to the state space and only depends on the previous point.
|
||
The properties of Markov processes have the following consequences:
|
||
\begin{align}
|
||
\sum_{x_{i+1}} p(\mathbf x_i \rightarrow \mathbf x_{i+1}) &= 1
|
||
\label{Markov1}
|
||
\end{align}
|
||
\begin{align}
|
||
\sum_{\mathbf x_i} \rho (\mathbf x_i) p(\mathbf x_i \rightarrow \mathbf x_{i+1}) &= \rho(\mathbf x_{i+1})
|
||
\label{Markov2}
|
||
\end{align}
|
||
the transition probability from configuration $i$ to configuration $j$ is $p(\mathbf x_i \rightarrow \mathbf x_j)$.
|
||
To ensure the validity of equation~\ref{Markov2}, the condition of microreversibility is sufficient that the movements in one direction are exactly compensated by the reverse movements.
|
||
This condition is verified by the detailed balance sheet equation:
|
||
\begin{align}
|
||
\rho(\mathbf x_i)p(\mathbf x_i \rightarrow \mathbf \mathbf x_j)=\rho(\mathbf x_j)p(\mathbf x_j \rightarrow \mathbf x_i)
|
||
\label{balanceEq}
|
||
\end{align}
|
||
The transition probability $p(\mathbf x_i \rightarrow \mathbf x_j)$ can be defined as the product of the probability of attempting a transition $\eta (\mathbf x_i \rightarrow \mathbf x_j)$ with the probability of accepting this same transition $\sigma (\mathbf x_i \rightarrow \mathbf x_j)$.
|
||
In the case of symmetrical movements, $\eta (\mathbf x_i \rightarrow \mathbf x_j)=\eta (\mathbf x_j \rightarrow \mathbf x_i)$ and the probability of accepting a relocation must satisfy the condition of the detailed assessment:
|
||
\begin{align}
|
||
\frac{\sigma (\mathbf x_i \rightarrow \mathbf x_{i+1})}{\sigma (\mathbf x_{i+1} \rightarrow \mathbf x_i)}=\frac{\rho (\mathbf x_{i+1})}{\rho(\mathbf x_i)}
|
||
\label{assess}
|
||
\end{align}
|
||
|
||
\subsection{Classical Molecular Dynamics}
|
||
|
||
%<<<<<<< HEAD
|
||
%MD is a powerful tool for analyzing the physical movements of atoms and molecules of many-body condensed matter systems. MD was originally developed following the earlier successes with Monte Carlo simulations. The first work about MD was published in 1957 by Alder \textit{et al.} in which the calculation results relate to the integration of classical motion equations for a system of hard spheres.\cite{Alder1957} Before long, the radiation damage at low and moderate energies were studied using MD in 1960 and the MD was simulated in liquid argon in 1964. \cite{Gibson1960, Rahman1964} MD experienced an extremely rapid development in the years that followed. MD simulation has been applied in chemistry, biochemistry, physics, biophysics, materials science, and branches of engineering, which is often coupled with experimental measurements to facilitate interpretation. MD has a strong predictive potential thus making it possible to motivate the implementation of new experiments. The diversity, broadness, and sophistication level of MD technique have been reported a lot.\cite{Allen2017, Cicccotti1987, Van1990, Berne1998, Tuckerman2000, Frenkel2001, Karplus2002, Rapaport2004, Tavan2005, Van2006} The field of the application of MD simulation is extremely wide, for instance, the study of structure, \cite{Parrinello1980, Ballone1988, Gutierrez1996, Karplus2002} thermodynamic, \cite{Briant1975, Postma1982,Honeycutt1987, Dang1997, Coveney2016} diffusion, \cite{Charati1998, Yeh2004, Braga2014, Pranami2015} viscosity, \cite{Mondello1997, Hess2002, Yeh2004, Zhang2015} spectroscopic. \cite{Aragon1976, Vuilleumier1999, Pavone2006, Brancato2006} Moreover, MD simulation is not only limited to the study of homogeneous systems but also allows the description of phase equilibria, the relaxation of metastable states and the dynamics of processes at interfaces. In addition, MD can also model chemical reactions in complex environments.
|
||
%
|
||
%MD simulations allow to do a real time evolution, one can then access certain properties explicitly dependent on time. The MD simulation can reach millisecond-scale.\cite{Shaw2009} The atoms or molecules are allowed to interact with each other in a fixed period of time, which gives a view of the dynamic evolution of the system. The evolution over time is calculated by solving number of Newton's motion equations for a system with interacting particles,
|
||
%For classical molecular dynamics, the forces are derived from a potential for empirical interaction between particles.
|
||
%A MD simulation needs the definition of a potential function, or a description of the terms by which the particles in the simulation will interact. This is often referred to as a force field in chemistry and biology, and as an interatomic potential in materials physics. Potentials can be defined at different levels of physical accuracy. The most commonly used in chemistry is based on molecular mechanics, which embodies a classical mechanics treatment of the interactions between particles that can reproduce conformational changes, however, it usually can not reproduce the chemical reactions. The reduction from a full quantum description to a classical potential involves two main approximations. The first one is the BO approximation (see section \ref{BO_approx}), which states that the dynamics of electrons are so fast that they can be considered to react instantaneously to the changes in the nuclei configuration, in other words the electronic wavefunction remains exactly in the variational ground state. As a consequence, the electrons and nuclei can be treated separately. The second one treats the nuclei (much heavier than electrons) as point particles that follow classical Newtonian dynamics. In a BOMD (Born-Oppenheimer molecular dynamics simulation, a full SCF calculation of the wavefunction is conducted at each timestep. The effect of the electrons is approximated as one PES, which usually represents the ground state.
|
||
%
|
||
%Classical MD is different from quantum dynamics approaches in which the temporal evolution of a system is described by the time-dependent Schr{\"o}dinger equation. For example, the multi-configuration time-dependent Hartree approach \cite{Beck2000, Wang2003} and discrete-variable representations, \cite{Wei1994, Light2000} which are particularly precise, but they can only treat a limited number of degrees of freedom.
|
||
%In this thesis, classical MD was used to perform simulations and we use the term MD to denote the classical MD.
|
||
|
||
%\textbf{Principles}
|
||
|
||
%In the classical formulation of MD, each particle or nucleus in the system is represented by a material point which interacts with all other particles via a potential defined by their positions. The principle of the dynamics of $N$ atoms is to determine the forces $\textbf F_\alpha$ acting on each of the particles in a given geometry and then calculate the accelerations and the velocities of the particles from these forces using the classical Newton's second law:
|
||
%=======
|
||
MD is a powerful tool for analyzing the physical movements of atoms and molecules of many-body systems.
|
||
MD was originally developed following the earlier successes of Monte Carlo simulations. The first work about
|
||
MD was published in 1957 by Alder \textit{et al.} which focused on the integration of classical equations of
|
||
motion for a system of hard spheres.\cite{Alder1957} Before long, radiation damage at low and moderate
|
||
energies were studied using MD in 1960 and MD was also applied to simulate liquid argon in 1964.\cite{Gibson1960, Rahman1964}
|
||
MD experienced an extremely rapid development in the years that followed. MD simulations have been applied
|
||
in chemistry, biochemistry, physics, biophysics, materials science, and branches of engineering, which is
|
||
often coupled with experimental measurements to facilitate interpretation. MD has a strong predictive
|
||
potential thus making it possible to motivate the implementation of new experiments. The diversity, broadness,
|
||
and sophistication level of MD techniques have been continuously reported.\cite{Allen2017, Cicccotti1987, Van1990, Berne1998, Tuckerman2000, Frenkel2001, Karplus2002, Rapaport2004, Tavan2005, Van2006} The range of applications of MD simulations is extremely wide, for instance,
|
||
the study of structure,\cite{Parrinello1980, Ballone1988, Gutierrez1996, Karplus2002}
|
||
thermodynamic,\cite{Briant1975, Postma1982,Honeycutt1987, Dang1997, Coveney2016}
|
||
diffusion,\cite{Charati1998, Yeh2004, Braga2014, Pranami2015}
|
||
viscosity,\cite{Mondello1997, Hess2002, Yeh2004, Zhang2015}
|
||
and spectroscopies.\cite{Aragon1976, Vuilleumier1999, Pavone2006, Brancato2006}
|
||
Moreover, MD simulations are not only limited to the study of homogeneous systems but also allow for
|
||
the description of phase equilibria, the relaxation of metastable states and the dynamics of processes
|
||
at interfaces. In addition, MD can also model chemical reactions in complex environments.
|
||
|
||
MD simulations allow to model real time evolution of particles, one can then access time-dependent properties. If the
|
||
time evolution is obtained by \textbf{integrating Newton's equations of motion} for a system of interacting particles, it is
|
||
referred to as \textbf{classical MD}.
|
||
A classical MD simulation needs the definition of a potential describing the interaction between particles in order to calculate
|
||
the PES. %This is often referred to as a force field in chemistry and biology, and as an interatomic potential in materials physics.
|
||
Potentials can be of different levels of accuracy as described in the previous section. The most commonly used potentials
|
||
in chemistry are \textbf{force fields}. In that case, one refers to \textbf{molecular mechanics}, which embodies a classical
|
||
mechanics treatment of the interactions between particles. In classical molecular dynamics, electrons and nuclei are not
|
||
distinguished and one refers only to particles. As already mentioned, the main drawback of force fields is that
|
||
they usually can not model chemical reactions. If the potential comes from a quantum chemical treatment of the
|
||
electrons only, the nuclei being treated as interacting point charge particles, on refers to \textbf{ab initio molecular dynamics}.
|
||
\textbf{Quantum dynamics} differs from classical molecular dynamics as the temporal evolution of a system is described by the
|
||
time-dependent Schr{\"o}dinger equation. For instance, the MCTDH approach\cite{Beck2000, Wang2003} and discrete-variable
|
||
representations,\cite{Wei1994, Light2000} are particularly accurate but are limited in terms of number of degrees of freedom.
|
||
The reduction from a full quantum description of all particles, electrons and nuclei, to a classical treatment involves two
|
||
main approximations. The first one is the BO approximation as described in section~\ref{BO_approx} which allows to
|
||
treat separately electrons and nuclei. The second one treats the nuclei (much heavier than electrons) as point charge particles
|
||
that follow classical Newtonian dynamics.
|
||
In this thesis, \textbf{classical molecular dynamics was used} to perform simulations and we use the term MD to denote
|
||
classical molecular dynamics only.
|
||
|
||
\textbf{Principles.}
|
||
In the classical formulation of MD, each particle or nucleus in the system is represented by a material point which
|
||
interacts with all other particles via a potential defined by their positions. The principle of the dynamics of $N$ atoms
|
||
is to determine the forces $\textbf F_\alpha$ acting on each of the particles in a given geometry and then calculate
|
||
the accelerations and the velocities of the particles from these forces using Newton's second law:
|
||
|
||
%>>>>>>> 92023a10c3aa8b7dc4ace43987c1d571fb99a738
|
||
\begin{align}
|
||
\mathbf F_\alpha=m_\alpha \mathbf a_\alpha, \ \ \ \alpha=1,2 ..., N
|
||
\label{Newton}
|
||
\end{align}
|
||
where $m_\alpha$ is the mass of the atom $\alpha$, $\textbf a_\alpha$ being its acceleration and $\textbf F_\alpha$ is the total force exerted
|
||
on $\alpha$. $\textbf F_\alpha$ is defined as the derivative of the potential energy $V$ of the system with respect to the corresponding
|
||
position ($\textbf R_\alpha$) of $\alpha$:
|
||
|
||
\begin{align}
|
||
\mathbf F_\alpha=- \frac{\partial V}{\partial \mathbf R_\alpha}=m_\alpha \frac{d^2 \mathbf R_\alpha}{dt^2}
|
||
\label{Force}
|
||
\end{align}
|
||
where $V$ depends on the positions of all atoms or particles. This leads to a system of $f \times N$ second-order differential equations
|
||
where $f$ is the dimension of space. In our case, the degree of freedom $f$ is equal to 3. According to the known initial positions of
|
||
particles, the potential energy can be obtained.
|
||
Then, a numerical resolution of the partial derivative equations provided in equation~\ref{Force} can be obtained using a suitable \textbf{integration algorithm}.
|
||
The integration algorithm gives access to the positions and velocities of atoms or particles and to the forces acting on these atoms
|
||
or particles over time. Here, the time is discretized in regular intervals and calculations are repeated at each time interval referred to as
|
||
\textbf{time step}.
|
||
%within the limit of the numerical integration step and the number of integrations carried out.
|
||
Many high order integration algorithms have been proposed depending on the desired accuracy:
|
||
Euler algorithm,\cite{Batina1991, Butcher2008} Verlet algorithms, predictor-corrector algorithm,\cite{Gear1966, Diethelm2002} and
|
||
Runge-Kutta algorithm.\cite{Runge1895, Kutta1901, Butcher1996, Butcher2015}
|
||
The Verlet algorithms include the Simple Verlet (SV),\cite{Verlet1967} the Leapfrog Verlet (LFV),\cite{Fraige2004} and the Velocity Verlet (VV).\cite{Swope1982}
|
||
The \textbf{Velocity Verlet algorithm} is the most widely used in many MD codes owing to its numerical stability and implementation simplicity.
|
||
Furthermore, movement constants are very well preserved over time. We briefly describe the basis of the Velocity Verlet algorithm below.
|
||
|
||
%The time symmetry in the Verlet algorithm reduces the level of local errors introduced into the integration by the discretization through removing all the odd-degree terms which refers to the three-degree terms in $\delta t$.
|
||
The local error is quantified by inserting the exact values $\textbf R_\alpha(t_{n-1})$, $\textbf R_\alpha(t_n)$, and $\textbf R_\alpha(t_{n+1})$ into the iteration and calculating the Taylor expansions at time $t=t_n$ of the position vector $\textbf R_\alpha(t\pm\delta t)$ in different time directions:
|
||
|
||
\begin{align}
|
||
\mathbf R_\alpha(t+\delta t) = \mathbf R_\alpha(t) + \mathbf v_\alpha(t)\delta t +\frac{\mathbf a_\alpha(t)\delta t^2}{2}+ \frac{\mathbf b_\alpha(t)\delta t^3}{6} + O(\delta t^4)
|
||
\label{DiscretisationError1}
|
||
\end{align}
|
||
|
||
\begin{align}
|
||
\mathbf R_\alpha(t-\delta t) = \mathbf R_\alpha(t) - \mathbf v_\alpha(t)\delta t +\frac{\mathbf a_\alpha(t)\delta t^2}{2} - \frac{\mathbf b_\alpha(t)\delta t^3}{6} + O(\delta t^4)
|
||
\label{DiscretisationError2}
|
||
\end{align}
|
||
where $\textbf{v}_\alpha(t)$ is the velocity of $\alpha$ and $\textbf{b}_\alpha(t)$ is the derivative of $\textbf{a}_\alpha(t)$ with respect to the time.
|
||
Through summing equations~\ref{DiscretisationError1} and \ref{DiscretisationError2}, we can get the Verlet integrator:
|
||
|
||
\begin{align}
|
||
\mathbf R_\alpha(t+\delta t) = 2\mathbf R_\alpha(t) - \mathbf R_\alpha(t-\delta t) + \mathbf a_\alpha(t)\delta t^2 + O(\delta t^4)
|
||
\label{Verlet}
|
||
\end{align}
|
||
|
||
We can notice that the first-order and third-order terms cancel out from the Taylor expansion, which makes the Verlet integrator more accurate
|
||
than the integration by a Taylor expansion only.
|
||
%It should also be noted that here the acceleration is calculated from the solution, but it is calculated at the central iteration point in the iteration process. %In the calculation of the global error, that is the distance between exact solution and approximation sequence, those two terms can not cancel out exactly which influences the global error order.
|
||
|
||
We can see in equation~\ref{Verlet} that the position propagation equation does not involve the velocities. They can be computed by the
|
||
following finite difference:
|
||
\begin{align}
|
||
\mathbf v_\alpha(t) = \frac{R_\alpha(t + \delta t) - R_\alpha(t - \delta t)}{2 \delta t} + O(\delta t^2)
|
||
\label{Velo}
|
||
\end{align}
|
||
|
||
This provides velocities at time $t$ and not at time $t+\delta t$, which means the velocity term is a step behind the position term.
|
||
The use of equation~\ref{Velo} has the advantage of low data storage, \textit{i.e.} less memory is required, but a problem emerges
|
||
in the calculation of the kinetic energy at time $ t+ \delta t$.
|
||
|
||
The VV algorithm is often applied to solve this problem as it allows velocities and positions to be computed simultaneously:
|
||
\begin{align}
|
||
\mathbf R_\alpha(t+\delta t) &= \mathbf R_\alpha(t)+ \mathbf v_\alpha(t)\delta t + \frac{\mathbf a_\alpha(t)\delta t^2}{2} \\
|
||
\mathbf v_\alpha(t+\delta t) &= \mathbf v_\alpha(t)+ \frac{\mathbf a_\alpha(t) + \mathbf a_\alpha(t+ \delta t)}{2} \delta t
|
||
\label{VV}
|
||
\end{align}
|
||
|
||
%The error in the VV is of the same order with the one in the basic Verlet. However, compared to the basic Verlet, the VV algorithm requires less memory because two vectors of position are tracked in basic Verlet while one position vector and one velocity vector are tracked in VV.
|
||
The standard implementation of the VV algorithm is a four steps scheme: firstly to calculate equation~\ref{VV1}, secondly to calculate equation~\ref{VV2}, thirdly to
|
||
derive $\mathbf a_\alpha(t+ \delta t)$ from the interaction potential with $\mathbf R_\alpha(t+ \delta t)$, finally to compute equation~\ref{VV3}.
|
||
|
||
\begin{align}
|
||
\mathbf v_\alpha(t+\frac{1}{2}\delta t) = \mathbf v_\alpha(t) + \frac{1}{2}\mathbf a_\alpha(t)\delta t
|
||
\label{VV1}
|
||
\end{align}
|
||
|
||
\begin{align}
|
||
\quad \qquad \mathbf R_\alpha(t+ \delta t) =\mathbf R_\alpha(t) + \mathbf v_\alpha(t + \frac{1}{2} \delta t) \delta t
|
||
\label{VV2}
|
||
\end{align}
|
||
|
||
\begin{align}
|
||
\qquad \qquad \qquad \mathbf v_\alpha(t+ \delta t) = \mathbf v_\alpha(t + \frac{1}{2} \delta t) + \frac{1}{2} \mathbf a_\alpha(t+ \delta t) \delta t
|
||
\label{VV3}
|
||
\end{align}
|
||
<%<<<<<< HEAD
|
||
%The VV algorithm needs to ensure two intrinsic properties of classical motion equations.
|
||
%One is temporal reversibility, the invariance of the trajectories at $t$ and $-t$. This symmetry leads to the independence of the dynamics of the microscopic system from the direction of time. The other one is the conservation of the Hamiltonian function over time. Because of the discretization of the trajectories, this conservation can not be assured. So a stable integration algorithm must impose this conservation for enough long time steps ($\delta t$) to allow the long simulation times. The VV algorithm is able to do this due to its sufficient numerical stability. At the end of each integration step, VV algorithm allows to access to $\mathbf R_\alpha (t+\delta t)$, $\mathbf v_\alpha (t+\delta t)$, $\mathbf F_\alpha (t+\delta t)$ directly.
|
||
|
||
%The force calculation in the MD simulations is an extremely significant part since it determines the precision level of the potential energy and the computation cost. Different methods were described in section \ref{electronicEnergy} to calculate the force. In the work of this thesis, SCC-DFTB method is applied to compute the force in the MD simulation.
|
||
|
||
%According to the conservation of energy, the natural set corresponding to this dynamic is the microcanonical set ($N, V, E$).\cite{Ray1981, Ray1999} $N$ is the number of particles. $V$ denotes the volume and $E$ is the energy of the system which is conserved. It is possible to extend the MD to other statistical sets by modifying the Hamiltonian of the system. For an example, it is the canonical set ($N, V, T$) when we add a thermostat to the system allowing to control the temperature $T$. \cite{Nose1984, Nose1984M, Hoover1985, Hunenberger2005}
|
||
Temperature is not a state variable in the simulations, an average kinetic temperature can be calculated from the average value of the kinetic energy.
|
||
|
||
%\subsection{\label{sec:PTMD} Parallel Tempering Molecular Dynamics}
|
||
|
||
%It is necessary to explore the energy surfaces as thoroughly as possible for the study the dynamical, thermodynamical, and structural properties of systems.
|
||
%For atomic and molecular system, PES usually presents a large number of stable configurations linked by energy barriers. However, it is not possible for MD simulations to overcome the energy barriers in a reasonable simulation time even with canonical set ($N, V, T$) the simulations are not ergodic. Actually, when a simulation explores the well of PES, it may have a high risk that this simulation will be blocked in this well at low temperature because the energy barriers are too high to be crossed. In this case, if $E_b$ refers to the energy barrier and $T$ is the temperature of the simulation, it can be consider that $E_b >> k_BT$. At an average temperature, the possibility of crossing the energy barriers during a simulation is increased. But it can not guarantee PES to be explored exhaustively. For high temperatures, it has a high probability to cross the barriers at high temperature, however, the bottoms of the well can not be explored comprehensively. Therefore, it is not possible to both pass the energy barrier and explore the bottom of the well exhaustively, only through verifying the $T$ in the MD simulation. Many methods have been proposed to solve this problem such as the enhanced sampling methods which are classified into two groups: biased methods and non-biased methods.
|
||
%For biased methods, the dynamics of the system are influenced by a external factor, usually a non-physical force which makes it possible to push the system outside the well even at low $T$, \cite{Torrie1977, Hansmann1993, Marchi1999, Bartels2000, Darve2001} for instance the Metadynamics. \cite{Laio2002, Iannuzzi2003, Barducci2011}
|
||
%For non-biased methods, the dynamics of the system are not modified directly. For instance, the simulated annealing, \cite{Kirkpatrick1983, Van1987} and multi-replica approaches such as the PTMD, which was used in this thesis.
|
||
%For simulated annealing method, a MD or Monte Carlo simulation is performed in the canonical ensemble by slowly decreasing the temperature of the thermostat in order to find a stable structure. The system is then heated again to be out from the basin that was just explored and then it is cooled again to find a new isomer. This mechanism should be done a lot of times in order to get the most stable structure.
|
||
%=======
|
||
%>>>>>>> 92023a10c3aa8b7dc4ace43987c1d571fb99a738
|
||
|
||
At the end of each integration step, the VV algorithm gives directly access to $\mathbf R_\alpha (t+\delta t)$, $\mathbf v_\alpha (t+\delta t)$,
|
||
and $\mathbf F_\alpha (t+\delta t)$.
|
||
The VV algorithm needs to ensure two intrinsic properties of the classical equations of motion. One is temporal reversibility,
|
||
the invariance of the trajectories at $t$ and $-t$. This symmetry leads to the independence of
|
||
he dynamics from the direction of time. The other property is the conservation of the Hamiltonian function over time.
|
||
Because of the discretization of trajectories, this conservation can not be insured. A stable integration algorithm must
|
||
impose this conservation for long enough time steps ($\delta t$) to allow for sufficiently long simulation times. The VV algorithm
|
||
is able to do this due to its sufficient numerical stability.
|
||
According to the conservation of energy, the natural ensemble corresponding to such dynamics is the microcanonical ensemble
|
||
($N, V, E$).\cite{Ray1981, Ray1999} $N$ is the number of particles. $V$ denotes the volume and $E$ is the energy of the system
|
||
which is conserved. It is possible to extend the MD to other statistical ensembles by modifying the Hamiltonian of the system.
|
||
For instance, in the canonical ensemble ($N, V, T$), a thermostat is added to the system allowing for the control of the temperature
|
||
$T$.\cite{Nose1984, Nose1984M, Hoover1985, Hunenberger2005} The temperature is not a state variable in the simulations,
|
||
an average kinetic temperature can be calculated from the average value of the kinetic energy.
|
||
|
||
The force, \textit{i.e.} acceleration, calculation in the MD simulations is an extremely important part as it determines both the accuracy
|
||
of the potential energy and the computational cost. Different methods were presented in section~\ref{electronicEnergy} to calculate the
|
||
forces. \textbf{In this thesis, the SCC-DFTB method has been applied to compute the forces in all the MD simulations. }
|
||
|
||
\subsection{Parallel-Tempering Molecular Dynamics} \label{sec:PTMD}
|
||
|
||
In a number of cases, it is necessary to explore the PES as thoroughly as possible for the study of dynamical,
|
||
thermodynamical, and structural properties of a given system. For atomic and molecular clusters, PES usually
|
||
presents a large number of stable configurations linked together by energy barriers. Unfortunately, MD simulations
|
||
cannot overcome these energy barriers in a reasonable simulation time even within the canonical ensemble ($N, V, T$).
|
||
This leads to non-ergodic simulations that cannot be used to extract meaningful statistical averages. Actually, when
|
||
a simulation explores the well of a PES, it may often be blocked in this well at low temperature because the energy
|
||
barriers are too high to be crossed. In this case, if $E_b$ refers to the energy barrier and $T$ is the temperature of
|
||
the simulation, one can consider that $E_b >> k_BT$ where $k_B$ is the Boltzmann constant. At intermediate temperature,
|
||
the possibility of crossing the energy barriers during a simulation increases, but this can not guarantee the PES to
|
||
be explored exhaustively. For high temperatures, one has a high probability to cross the energy barriers whereas the
|
||
bottoms of the wells can not be explored comprehensively. Therefore, it is not possible to both cross the energy barriers
|
||
and thoroughly explore the bottom of the wells using a unique MD simulation at a given temperature.
|
||
Many methods have been proposed to solve this question and are referred to as \textbf{enhanced sampling methods}
|
||
They are classified into two groups: \textbf{biased methods} and \textbf{non-biased methods}.
|
||
In \textbf{biased methods}, the dynamics of the system is influenced by a external factor, usually a non-physical
|
||
orce, which makes it possible to push the system outside of the wells even at low $T$.\cite{Torrie1977, Hansmann1993, Marchi1999, Bartels2000, Darve2001}
|
||
For instance, Metadynamics is a biased method.\cite{Laio2002, Iannuzzi2003, Barducci2011}
|
||
In \textbf{non-biased methods}, the dynamics of the system is not modified directly. Examples are simulated annealing,\cite{Kirkpatrick1983, Van1987}
|
||
and multi-replica approaches such as the \textbf{parallel-tempering molecular dynamics} approach, which has been used in this thesis.
|
||
%For simulated annealing method, a MD or Monte Carlo simulation is performed in the canonical ensemble by slowly decreasing the temperature of the thermostat in order to find a stable structure. The system is then heated again to be out from the basin that was just explored and then it is cooled again to find a new isomer. This mechanism should be done a lot of times in order to get the most stable structure.
|
||
|
||
The replica exchange approach also termed parallel-tempering was originally devised by Swendsen $et~ al.$ in 1986 \cite{Swendsen1986}
|
||
then extended by Geyer and coauthor in 1991 \cite{Geyer1991} and was further developed by Hukushima and Nemoto,\cite{Hukushima1996}
|
||
Falcioni and Deem,\cite{Falcioni1999} Earl and coworker.\cite{Earl2005} Sugita and Okamoto formulated a MD version of parallel tempering
|
||
to enhance conformational sampling.\cite{Sugita1999} PTMD is a method which aims at enhancing the ergodicity of MD simulations.
|
||
The principle of PTMD is shown in Figure~\ref{ptmd_s}.
|
||
|
||
\begin{figure}[h]
|
||
\begin{center}
|
||
\includegraphics[width=12cm]{PTMD_schema.png}
|
||
\caption{Schematics of the PTMD algorithm in its synchronous version. Replicas of
|
||
the same system, numbered from $C_1$ to $C_N$, are simulated subject to different temperatures (from T$_1$ to T$_N$).
|
||
Once $X$ MD steps (straight solid arrows) have been performed by each replica, configuration
|
||
exchanges are attempted between neighbouring simulations according to the Metropolis criterion.
|
||
Some of them undergo successful swapping (solid curved arrows) while other not (dashed curved arrows).
|
||
MD simulations then proceed for $X$ additional MD steps before new attempts of exchange.} \label{ptmd_s}
|
||
\end{center}
|
||
\end{figure}
|
||
|
||
$N$ replicas ($C_i, i = 1, 2, ..., N$) of the same system are simulated in parallel each one at a given temperature $T_i,$
|
||
(i = 1, 2, ..., N) in the canonical ensemble.
|
||
%The configurations at high $T$ can be available to be simulated at low $T$ and the configurations at low $T$ also have access to the simulations at high $T$.
|
||
%This made it possible to model configurations at both high and low temperatures through PTMD simulations.
|
||
The time evolution of each replica is independent with each other but exchanges of configurations between adjacent
|
||
replicas $C_i$ and $C_j$, where $T_i < T_j$ and $i = j - 1$ are permitted at regular time intervals.
|
||
The choice of the extreme temperatures $T_1$ and $T_N$ is very important for the algorithm to be optimal.
|
||
The lowest temperature ($T_1$) should be the one at which usual simulations are blocked and the highest
|
||
temperature ($T_N$) should be chosen so that all significant energy barriers can be overcome during the
|
||
simulation. Moreover, the temperatures between $T_1$ and $T_N$ must be chosen to lead to sufficient overlap
|
||
between the density of states of the adjacent replicas. Indeed, if this overlap is too small, the probability of
|
||
exchange is very low, which makes the PTMD simulations inefficient and leads to a bad exploration of the PES.
|
||
In contrast, if the overlap is too large, a large amount of redundant information will be produced, which will
|
||
cost unnecessary computational resources. Configurations between two neighbouring replicas at different
|
||
$T$ are exchanged based on the Metropolis–Hastings criterion with probability:
|
||
\begin{align}
|
||
\rho (C_i \Leftrightarrow C_j)=\mathrm{min} \left(1, e^{(V_i-V_j) \left(\frac{1}{kT_i} - \frac{1}{kT_j} \right)} \right) = \mathrm{min} \left(1, e^m \right)
|
||
\label{Metropolis}
|
||
\end{align}
|
||
where $V_i$ and $V_j$ are the potential energies of replicas $i$ and $j$, respectively. $T_i$ and $T_j$ are the temperatures
|
||
of replicas $i$ and $j$, respectively.
|
||
If the energy of replica configuration $C_j$ at high temperature $T_j$ is lower than the energy of replica configuration
|
||
$C_i$ at low temperature $T_i$, which means the exponent ($m$) of $e$ is positive, then the exchange is allowed.
|
||
If $m$ is negative, the exchange between neighbouring replicas is only allowed when $e^m$ is greater than $\omega$
|
||
(a random value between 0 and 1). To accelerate the equilibration of the system after the exchange, the velocities of
|
||
all particles can be renormalized as follows:
|
||
\begin{align}
|
||
\mathbf v_\alpha^\mathrm{new}=\left( \frac{T_\mathrm{new}}{T_\mathrm{old}}\right)^{\frac{1}{2}} \mathbf v_\alpha^\mathrm{old}
|
||
\label{newVelo}
|
||
\end{align}
|
||
|
||
%Generally, the thermodynamical properties can not be well calculated with stand MD simulations while they can be computed
|
||
%with fairly good accuracy through PTMD simulations. It should be noted that the heat capacity can be obtained through PTMD
|
||
%simulations but it is not available through simulated annealing or basin hopping.
|
||
%The heat capacity curves of aggregates have been certificated that they can provide important information to the behavior (such as phase transition) of these aggregates as a function of temperature. \cite{Dougherty1998, Douady2009, Schmidt2012, Korchagina2017, Rapacioli2018}
|
||
|
||
|
||
\subsection{Global Optimization}
|
||
|
||
%<<<<<<< HEAD
|
||
%It usually includes a large number of stationary points in the PES, so it is not easy to find the global minimum (the lowest minimum). Only local optimization methods do not make it possible to cross the potential barriers between local minima to obtain the global minimum. Therefore, a global optimization method such as MD and Monte Carlo simulations is needed to perform a larger exploration of the energy surface to get the lowest minimum.
|
||
%=======
|
||
%>>>>>>> 92023a10c3aa8b7dc4ace43987c1d571fb99a738
|
||
|
||
\textbf{Global Optimization} refers to the determination of the lowest energy point on a PES, \textit{i.e.} the global minimum. As this latter usually
|
||
includes a large number of stationary points, it is not straightforward to find the global minimum. Local optimization methods do not
|
||
make it possible to cross the energy barriers between local minima. Therefore, a global optimization scheme such as MD or Monte Carlo
|
||
simulations is needed to perform a more exhaustive exploration of the PES to get to the lowest energy minimum.
|
||
|
||
There exits a vast amount of methods to perform global optimization and each one has its strength and weaknesses.
|
||
For instance, the Basin-Hopping method is a particular useful global optimization technique in high-dimensional landscapes
|
||
that iterates by performing a random perturbation of coordinates, making a local optimization, and rejecting or accepting
|
||
new coordinates based on a minimized function value.\cite{Wales1997, Wales1999}
|
||
Genetic algorithms are also among the most used methods to find a global minimum.\cite{Hartke1993, Unger1993, Sivanandam2008, Toledo2014}
|
||
A genetic algorithm is inspired by the process of natural selection. Genetic algorithms are usually applied to generate high-quality solutions
|
||
of optimization. The ergodicity problem appears in all of these global optimization methods. In principle, one can only be sure of having
|
||
found the real global minimum after an infinite number of iterations.
|
||
|
||
In order to avoid ergodicity problems, it is interesting to combine global and local optimization methods. A very popular combination
|
||
is the simulated annealing method combined with local optimizations. The PES of molecular aggregates display many degrees of freedom
|
||
and contains a large number of low-energy isomers. The PTMD algorithm coupled with a great number of local optimizations is a good
|
||
choice to search for low-energy structures in this kind of system. Local optimizations are performed many times from initial conditions
|
||
structures which are extracted from all PTMD trajectories, whether it be low or high temperature, in order to maximize sampling. \textbf{This
|
||
approach, in combination with SCC-DFTB, has been conducted along this thesis to perform global Optimization.}
|
||
|
||
|
||
|
||
|
||
|
||
%: ----------------------- HELP: special characters
|
||
% above you can see how special characters are coded; e.g. $\alpha$
|
||
% below are the most frequently used codes:
|
||
%$\alpha$ $\beta$ $\gamma$ $\delta$
|
||
|
||
%$^{chars to be superscripted}$ OR $^x$ (for a single character)
|
||
%$_{chars to be suberscripted}$ OR $_x$
|
||
|
||
%> $>$ greater, < $<$ less
|
||
%≥ $\ge$ greater than or equal, ≤ $\ge$ lesser than or equal
|
||
%~ $\sim$ similar to
|
||
|
||
%$^{\circ}$C ° as in degree C
|
||
%± \pm plus/minus sign
|
||
|
||
%$\AA$ produces Å (Angstrom)
|
||
|
||
|
||
%: ----------------------- HELP: references
|
||
% References can be links to figures, tables, sections, or references.
|
||
% For figures, tables, and text you define the target of the link with \label{XYZ}. Then you call cross-link with the command \ref{XYZ}, as above
|
||
% Citations are bound in a very similar way with \cite{XYZ}. You store your references in a BibTex file with a programme like BibDesk.
|
||
|
||
|
||
%%\figuremacro{largepotato}{A common glucose polymers}{The figure shows starch granules in potato cells, taken from
|
||
%%\href{http://molecularexpressions.com/micro/gallery/burgersnfries/burgersnfries4.html}{Molecular Expressions}.}
|
||
|
||
%: ----------------------- HELP: adding figures with macros
|
||
% This template provides a very convenient way to add figures with minimal code.
|
||
% \figuremacro{1}{2}{3}{4} calls up a series of commands formating your image.
|
||
% 1 = name of the file without extension; PNG, JPEG is ok; GIF doesn't work
|
||
% 2 = title of the figure AND the name of the label for cross-linking
|
||
% 3 = caption text for the figure
|
||
|
||
%: ----------------------- HELP: www links
|
||
% You can also see above how, www links are placed
|
||
% \href{http://www.something.net}{link text}
|
||
|
||
%%\figuremacroW{largepotato}{Title}{Caption}{0.8}
|
||
% variation of the above macro with a width setting
|
||
% \figuremacroW{1}{2}{3}{4}
|
||
% 1-3 as above
|
||
% 4 = size relative to text width which is 1; use this to reduce figures
|
||
|
||
%Insulin stimulates the following processes:
|
||
|
||
%%\begin{itemize}
|
||
%%\item muscle and fat cells remove glucose from the blood,
|
||
%%\item cells breakdown glucose via glycolysis and the citrate cycle, storing its energy in the form of ATP,
|
||
%%%%\item liver and muscle store glucose as glycogen as a short-term energy reserve,
|
||
%%\item adipose tissue stores glucose as fat for long-term energy reserve, and
|
||
%%\item cells use glucose for protein synthesis.
|
||
%%\end{itemize}
|
||
|
||
%: ----------------------- HELP: lists
|
||
% This is how you generate lists in LaTeX.
|
||
% If you replace {itemize} by {enumerate} you get a numbered list.
|
||
|
||
|
||
%: ----------------------- HELP: tables
|
||
% Directly coding tables in latex is tiresome. See below.
|
||
% I would recommend using a converter macro that allows you to make the table in Excel and convert them into latex code which you can then paste into your doc.
|
||
% This is the link: http://www.softpedia.com/get/Office-tools/Other-Office-Tools/Excel2Latex.shtml
|
||
% It's a Excel template file containing a macro for the conversion.
|
||
|
||
%%\begin{table}[htdp]
|
||
%%\centering
|
||
%%\begin{tabular}{ccc} % ccc means 3 columns, all centered; alternatives are l, r
|
||
|
||
%%{\bf Gene} & {\bf GeneID} & {\bf Length} \\
|
||
% & denotes the end of a cell/column, \\ changes to next table row
|
||
%%\hline % draws a line under the column headers
|
||
|
||
%% human latexin & 1234 & 14.9 kbps \\
|
||
%% mouse latexin & 2345 & 10.1 kbps \\
|
||
%% rat latexin & 3456 & 9.6 kbps \\
|
||
% Watch out. Every line must have 3 columns = 2x &.
|
||
% Otherwise you will get an error.
|
||
|
||
%%\end{tabular}
|
||
%%\caption[title of table]{\textbf{title of table} - Overview of latexin genes.}
|
||
% You only need to write the title twice if you don't want it to appear in bold in the list of tables.
|
||
%%\label{latexin_genes} % label for cross-links with
|
||
%%\ref{latexin_genes}
|
||
%%\end{table}
|
||
% ---------------------------------------------------------------------- |