methods.tex

%methods

\acresetall

\part{Methods}
\label{sec:methods}

\chapter{Program Package \textsc{Spheres}}
\label{sec:theprogramspheres}

For the projects outlined in
chapters~\ref{sec:fromstickyhardspheretoLJtypeclusters} and
\ref{sec:thegregorynewtonclusters} a program to optimise 3-dimensional
coordinates of cluster structures with respect to a two-body potential was
required to be created. The optimisation routine needed to be flexible in the
way that it would be easy to implement different two-body potentials like
\acf{LJ} or \acf{eLJ}. Furthermore, the program needed to be able to analyse the
results regarding structure, energy and the matrix of second derivatives. The
resulting program was written in \Cpp with standard library version 11 and was
tested to compile in a Linux environment with the \texttt{clang++} compiler
version 6.0. In the following sections, main features are explained and more
important functionalities are outlined in more detail.


\section{Structural Optimisation and Analysis}
\label{sec:generalstructure}

The main functionality of the program is implemented in three different
executables. Each carries out a different task, i.e. optimising given input
structures and output the results, removing duplicate structures from a set of
input structures and finding differences of two sets of input structures. The
programs are set up in such a way that their outputs can be used as input files
for the other program parts. This allows the programs to be used sequentially
while retaining the flexibility of using each executable separately. This
chain-like execution scheme is depicted in
figure~\ref{fig:OptimisationAnalyseScheme}. 
%
\begin{figure}[htb]\centering
    \begin{tikzpicture}[node distance = 3cm, scale=0.8, auto, transform shape]%, every node/.style={scale=0.8}]
    \node [cloud] (Input1) {Input structures};
    \node [block, right of=Input1] (Optimise) {Optimise};
    \node [decision, right of=Optimise] (Hessian) {Negative values in Hessian?};
    \node [block, below of=Optimise] (reopt) {Displace geometry};
    \node [cloud, right=1.5cm of Hessian] (print) {Optimised structures};
    \node [block, below of=print] (Analyse) {Analyse $\{r_{ij}\}$};
    \node [cloud, below of=Analyse] (Analysed) {Duplicates removed};
    \node [block, left of=Analysed] (Match) {Match sets};
    \node [cloud, left of=Match] (Input2) {Input structures};
    \node [cloud, below of=Match] (missing) {Missing clusters};

    \path [line] (Input1) -- (Optimise);
    \path [line] (Optimise) -- (Hessian);
    \path [line] (Hessian) |- node {yes} (reopt);
    \path [line] (reopt) -- (Optimise);
    \path [line] (Hessian) -- node {no} (print);
    \path [line] (print) -- (Analyse);
    \path [line] (Analyse) -- (Analysed);
    \path [line] (Analysed) -- (Match);
    \path [line] (Input2) -- (Match);
    \path [line] (Match) -- (missing);

    \node (OptProg) [transform shape=false,draw, color=blue, dashed, thick, %
        fit=(Optimise) (Hessian) (reopt)] {};
    %\node [above=0.1cm of OptProg] {\texttt{Optimize}};

    \node (InputBox1) [draw,color=gray,thick,transform shape=false, fit=(Input1)] {}; 
    \node (InputBox2) [draw,color=gray,thick,transform shape=false, fit=(Input2)] {}; 
    \node (OutputBox1) [draw,color=gray,thick,dashed,transform shape=false, fit=(print)] {}; 
    \node (OutputBox2) [draw,color=gray,thick,dashed,transform shape=false, fit=(Analysed)] {}; 
    \node (OutputBox3) [draw,color=gray,thick,dashed,transform shape=false, fit=(missing)] {}; 
    \node [above=0.1cm of InputBox1] {Input 1};
    \node [below=0.1cm of InputBox2] {Input 2};
    \node [above=0.1cm of OutputBox1] {Output 1};
    \node [below=0.1cm of OutputBox2] {Output 2};
    \node [below=0.1cm of OutputBox3] {Output 3};

\end{tikzpicture}
    \caption{Schematic representation of the optimisation, analysis and
    matching procedure. Red circles: program input or output, blue ellipses:
    program executions.}
    \label{fig:OptimisationAnalyseScheme}
\end{figure}
%
A set of input structures (Input 1) is provided to the first executable program
\texttt{optimize}, which optimises the structures and generates Output 1, which
is a list of optimised coordinates in the same format as the input. This can be
used as input for the second executable program \texttt{analyze}, which uses
the distance matrix to identify duplicate structures and outputs a list of
optimised structures with no duplicates (Output 2). In combination with a
secondary set of structures (Input 2) the third executable program
\texttt{match} can be used to compare both sets of structures and output
coordinates of structures that are missing from either set 1 or 2 (Output 3). 

\subsection{Optimisation of Input Structures}
\label{sec:optimisationofinputstructures}

The optimisation of input structures with a chosen potential can be carried out
with the program \texttt{optimize}. The coordinates of the input structures have
to be provided in a single file where each line contains three numbers
representing the position of one sphere. Multiple structures can be included by
separating the list of coordinates by single blank lines. The program
automatically moves the input structures to their respective centres of mass,
and rotates them onto the principal axis system. To speed up the optimisation
the environment variable \texttt{OMP\_NUM\_THREADS} can be set to a value
greater than one to enable \ac{OMP} parallelisation. Each optimisation is
carried out in a separate thread, while the number of simultaneously running
threads is equal to \texttt{OMP\_NUM\_THREADS}.

Various parameters of the program run can be controlled using a
\texttt{settings} file in the working directory, which is parsed with help of
the library
\textit{libconfig}.\autocite{Lindner_libconfiglibraryprocessing_2018} For
example, the optimisation can be controlled using the \texttt{opt} group. In
this group, the optimisation model can be chosen with the \texttt{name} setting,
which currently can be set to either \texttt{BFGS} for the \acf{BFGS} algorithm
or \texttt{CG} for the conjugate gradient method. The optimiser uses the machine
learning library dlib\autocite{King_DlibmlMachineLearning_2009} as a back-end,
which allows for the implementation of additional optimisation models with
relative ease. The energy termination criterion is defined using the
\texttt{dforce} setting and the maximum number of steps can be set with the
\texttt{nsteps} setting. An example is given in the box below.
%\begin{figure}[ht]
%
\begin{Verbatim}[frame=single,label=settings file - opt tag]
opt: {
    name = "BFGS";
    dforce = 1e-10;
    nsteps = 100;
};
\end{Verbatim}
%\end{figure}
%
A potential model has to be chosen for the optimisation procedure. Custom models
can be added easily, which will be explained later. The current implementation
allows the selection of three pair potentials, \ac{LJ}, \ac{eLJ} and a \ac{LJ}
potential with a cut-off range. For the \ac{LJ} potential four parameters have
to be given to the \texttt{potential} group as shown in the box below.
%
%\begin{figure}[htb]
\begin{Verbatim}[frame=single,label=settings file - potential tag]
potential: {
    name = "LJ";
    epsilon = 1.0;
    rm =  1.0;
    exp1 = 12.0;
    exp2 = 6.0;
};
\end{Verbatim}
%\end{figure}
%
The \texttt{name} setting enables the \ac{LJ} potential and the four parameters
\texttt{epsilon}, \texttt{rm}, \texttt{exp1} and \texttt{exp2} define the
potential. The values shown in the example are the default values, which the
program will fall back to if no values are provided by the user. The \ac{eLJ}
potential can be chosen by setting the \texttt{name} setting to \texttt{ELJ} and
a range cut-off value can be chosen with the \texttt{RangeLJ} potential. For the
\ac{eLJ} potential the $c_n$ coefficients have to be provided in a separate file
located in the working directory and named \texttt{ext}, where each line
contains two numbers. The first integer represents $n$ and the second floating
point number defines the corresponding coefficient.

After each structural optimisation the program checks if the result is a true
minimum by calculating the Hessian matrix eigenvalues. If this check fails the
eigenvectors of the Hessian matrix are calculated and the non-minimum structure
is displaced in both possible directions according to the eigenvector belonging
to the first negative eigenvalue in the Hessian matrix. The algorithm tries to
re-optimise until there are no negative eigenvalues in the Hessian matrix or
until a maximum number of five re-optimisations is reached.

The resulting structures are printed in the file \texttt{coord} in the same
format as the input. To analyse the optimisation procedure the additional files
\texttt{opt} and \texttt{reopt} are created, which contain the intermediate
coordinates of the optimisation and re-optimisation procedures for each individual
input structure.

\subsection{Removing Duplicate Structures}
\label{sec:analysingresults}

Duplicated geometries in a set of input structures can be identified using the
second program \texttt{analyze}. The input has to be provided in the same format
as for the \texttt{optimize} program. The read-in procedure is equivalent to the
\texttt{optimize} program. The program uses two methods to identify unique
structures, one of which uses the energy of the cluster as a criterion,
therefore the \texttt{potential} group has to be set in the \texttt{settings}
file.

The first method uses four values to uniquely identify a structure. Those are
the values for energy and the three eigenvalues of the moment of inertia tensor.
This sorting procedure uses a \verb|map| container from the \Cpp standard
library, which is an implementation of a binary search tree. The \verb|map|
stores key-value pairs, which are guaranteed to be stored in order with respect
to the key.\autocite{Stroustrup_programminglanguage_2000} In this specific case,
the ordering is implemented as a custom comparator function that sorts by energy
first, then smallest moment of inertia tensor eigenvalue followed by the two
larger moment of inertia tensor eigenvalues. The comparator function returns
\texttt{true} if the first argument is considered smaller than the second
argument. This procedure is shown as pseudocode in algorithm~\ref{algo:sortEI}. 
%
\begin{algorithm}[htp]
    \caption{Comparator function for sorting by energy $E$ and moment of inertia eigenvalues $I_1\leq I_2\leq I_3$.}
    \label{algo:sortEI}
    \begin{algorithmic}[1]
        \Procedure{Compare}{structure $a$, structure $b$}
            \If{$E^a < E^b$} \State \textbf{return} true
            \ElsIf{$E^a = E^b$ \textbf{and} $I_1^a < I_1^b$} \State \textbf{return} true
            \ElsIf{$E^a = E^b$ \textbf{and} $I_1^a = I_1^b$ \textbf{and} $I_2^a < I_2^b$} \State \textbf{return} true
            \ElsIf{$E^a = E^b$ \textbf{and} $I_1^a = I_1^b$ \textbf{and} $I_2^a = I_2^b$ \textbf{and} $I_3^a < I_3^b$} \State \textbf{return} true
            \Else \State \textbf{return} false
            \EndIf
        \EndProcedure
    \end{algorithmic}
\end{algorithm}
%
The key is mapped to a value, which is simply the number of structures.
The advantage of this method is its great scalability with respect to the number
of input structures as values can be retrieved quickly based on the key.

The second method uses the \ac{EDM} as the differentiation criterion. The
\ac{EDM} is the matrix of all inter-particle distances $d_{ij}$ where
each entry is defined as the Euclidean norm $||\cdot||$ between two spheres.
%
\begin{align}
    d_{ij}=||\mathbf{x}_i-\mathbf{x}_j||^2
\end{align}
%
Each unique embedding of the cluster in space can be represented by an \ac{EDM},
however information about the absolute position, orientation and chirality is
not contained in this representation. That means rigid transformations of
clusters (translations, rotations, reflections) don't affect the \ac{EDM} as
they don't change fixed distances between points in
space.\autocite{Dokmanic_EuclideanDistanceMatrices_2015} If at least one
inter-particle distance is different, the structures are said to be not equal.
The algorithm is implemented in such a way that a structure's distance matrix is
compared to the distance matrix of all other already sorted structures and is
added to the matching group if they are equal up to a set threshold. If not, the
structure is added to the array as a new group. If the optimisation results in
many unique structures, this method of sorting becomes slow, as each trial
structure has to be compared to all other already sorted structure groups.

\subsection{Matching Structures}
\label{sec:matchingstructures}

To compare the results from the optimisations procedures to previously
published sets of clusters the program \verb|match| can be used. It takes two
files as input that each contain a set of structures of equal or different size
and compares them based on the \ac{EDM}. The number of atoms in each set must
be equal, otherwise the program will be terminated. If both sets are found to
be identical no output files are created. In case there are unmatched
structures they will be printed in \verb|xyz| format in the \verb|output|
directory for further analysis.

\section{Graph-Theoretical Analysis}
\label{sec:GraphTheoreticalAnalysis}

Graph theoretical analysis of cluster structures can be done with
\verb|ico-subgraph|, which is a very specialised program designed to compare
contact graphs of \ac{GN} clusters to the icosahedral graph. The input format
for the structures remains unchanged.\footnote{This part of the program was used
in chapter~\ref{sec:thegregorynewtonclusters}. For a definition of contact
graphs and \ac{GN} clusters refer to the introduction in that part of the
thesis.}

All graph objects are handled by the \textit{Boost Graph
Library}\autocite{Siek_BoostGraphLibrary_2002} and are implemented as undirected
graphs as the direction of the connectivity in clusters is not meaningful. The
graphs are automatically generated from three-dimensional coordinate input given
by an object of type \verb|structure|. Two types of graphs can be generated from
these: (1) the graph of the complete structure and (2) the contact graph
containing all spheres but the central one (if it exists). The latter is the
more important one as it was used to carry out the analyses in
chapter~\ref{sec:thegregorynewtonclusters}. The decision whether two spheres are
connected is based on two parameters: the equilibrium distance $r_e$ and a
threshold value $\epsilon$. The default values for these are $r_e=1$ and
$\epsilon=10^{-10}$, respectively.

\begin{align}
    d_{ij}-r_e < \epsilon
\end{align}

The icosahedral graph is compared to the input structures via their graphs and
the \textit{VF2} graph matching
algorithm.\autocite{Cordella_ImprovedAlgorithmMatching_2001,Cordella_SubGraphIsomorphism_2004}
The algorithm finds all mappings of the vertices of the icosahedral graph to the
graph of the input structure. As the icosahedron represents the complete planar
graph for 12 vertices, every graph that is a subgraph of the icosahedral graph
can be represented by the number and type of edges removed from the icosahedral
graph. For this application the mapping was chosen based on the \ac{RMS} value
of the distances between the spheres corresponding to removed edges. From all
the possible mappings of the investigated graphs to the icosahedral graph the
one with the lowest \ac{RMS} value was chosen.

The graphs of the input structures are analysed with respect to their vertex
and face degrees. Vertex degrees are calculated directly by the library and can
be accessed with the \verb|degree| function. For face degrees the
\verb|planar_face_traversal| function has to be invoked. This algorithm
iterates over all faces in the planar embedding counting the number of vertices
constructing each face. The graphs of the input structures are then sorted
based on the calculated degree values, starting with the vertex degree in
descending order and followed by the face degree in ascending order. If two or
more graphs have the same amount of face and vertex degrees they are grouped
together. The sorted graphs are printed to the standard output in the same form
as table~\ref{tab:verticesandfaces} (page~\pageref{tab:verticesandfaces}).
Additionally, the investigated graphs are printed in terms of removed edges
from the icosahedral graph as shown in table~\ref{tab:icosubgraphs}
(page~\pageref{tab:icosubgraphs}).


\section{Additional Functionalities}
\label{sec:AdditionalFunctions}

Besides the main parts of the program, which have been described on the previous
pages, a few script-like executable programs are provided. These were used to
calculate various different properties of the investigated clusters.

\paragraph{Analyse Bond Lengths}%app-bondvariance app-shortestbond (app-N14?)
The bond variance in optimised structures can be calculated with
\texttt{app-bondvariance}. The bond variance is simply defined as the difference
between the shortest and longest bond of a cluster structure. In clusters
optimised by soft potentials a bond is not as well defined as for \ac{SHS}
clusters. Therefore, it has to be determined with respect to a threshold value
and the variance of the bond lengths can not be larger than this threshold. This
application has been used in
chapter~\ref{sec:612LennardJonesClustersfromBasinHopping} to calculate the bond
variance of the optimised $(6,12)$-\ac{LJ} structures. A more specialised
version called \texttt{app-shortestbond} is also provided that can be used to
find the cluster with the shortest bond distance.

\paragraph{Sort Structures by $N_c$}%app-GN and app-Nc

Analysing the total contact numbers $N_c$ or specific kissing numbers can be
done with the programs \texttt{app-Nc} and \texttt{app-GN}. The latter looks
for clusters with a central atom that has exactly 12 spheres arranged around
it, so called Gregory-Newton clusters.

\section{Implementations in Detail}
\label{sec:ImplementationsinDetail}

In the following sections the basic implementation of cluster structures and
two-body potentials will be explained in more detail.

\subsection{Treating Cluster Structures}
\label{sec:thestructureclass}

The main purpose of program \textsc{Spheres} was to optimise an input set of
structures with given Cartesian coordinates and analyse the results of the
optimisation based on properties of the resulting structures. The handling of
those structures and respective properties was therefore crucial to the
functioning of the program.

The Cartesian coordinates are read from a file with blank lines separating
different starting geometries. For each of these individual sets of coordinates
a \verb|structure| object is created. The \verb|structure| class, however, has
much more functionality than storing 3D coordinates. In fact, it serves as a
complex data type storing properties besides coordinates as well as member
functions to calculate those properties. They are stored in data members defined
in private fields, such that they can only be manipulated by functions owned by
the \verb|structure| class. The \verb|structure| class is designed to store
values for:
\begin{itemize}
    \item Cartesian coordinates
    \item energy (depending on the chosen potential), 
    \item an integer number for labelling, 
    \item moment of inertia, 
    \item the \ac{EDM} and adjacency matrix representations and
    \item an \verb|undirectedGraph| that contains the connectivity information.
\end{itemize}

The Cartesian coordinates of each individual sphere are stored in an object of
type \texttt{coord3d}, which is a modified version of the \verb|coord3d|
implementation used in program
\textsc{Fullerene}.\autocite{Schwerdtfeger_Programfullerenesoftware_2013} 

The moment of inertia tensor determines the torque needed to accelerate a rigid
body to spin around a rotational axis through the origin of the coordinate
system. It is therefore analogous to mass in case of linear, translational
acceleration. For cluster structures the tensor is equal to the sum over the
moments of inertia of all constituent particles. As the clusters investigated in
this thesis are only made up of one type of particle the mass term was set to
unity. The inertia tensor $\mathbf{I}$ can then be calculated via the equations
below.
%
\begin{gather}
    \begin{gathered}
    \mathbf{I}=
    \begin{pmatrix}
        I_{xx} & I_{xy} & I_{xz}\\
        I_{yx} & I_{yy} & I_{yz}\\
        I_{zx} & I_{zy} & I_{zz}
    \end{pmatrix}\\
    I_{xx}=\sum_i(y_i^2+z_i^2) \quad I_{yy}=\sum_i(x_i^2+z_i^2) \quad I_{xx}=\sum_i(y_i^2+x_i^2)\\
    I_{xy} = I_{yx} = -\sum_ix_iy_i \quad I_{xz} = I_{zx} = -\sum_ix_iz_i \quad I_{yz} = I_{zy} = -\sum_iy_iz_i
    \end{gathered}
\end{gather}
%
Here, $x_i$, $y_i$ and $z_i$ denote the respective coordinate of sphere $i$.
Diagonalising the inertia tensor yields a set of eigenvalues and eigenvectors,
with the latter representing the principal axis system. Upon creation of a
\verb|structure| object the coordinates are transformed, such that the
coordinate origin lies at the centre of mass and the structures principal axis
are aligned with the basis vectors of the Cartesian coordinate system.

The class is designed to ensure that when any of the particle coordinates
change, all properties are recalculated such that there is never a mismatch
between the properties and the coordinates they refer to.


\subsection{Treating Two-Body Potentials}
\label{sec:thepairpotentialclass}

For the geometry optimisation of the cluster structures, methods to calculate the
energy and gradient need to be provided to the optimisation library. As shown in
section~\ref{sec:PracticalImplementationForPotentialsDependingOnPairDistances}
the gradient can be expressed as in equation~\eqref{eqn:gradient}.
%
\begin{align}
    \pdv{E(\mathbf{X})}{x_m}=\sum_{j>i}^N\pdv{\varepsilon(r_{ij})}{r_{ij}}\ \pdv{r_{ij}}{\mathbf{r}_{ij}}\ \pdv{\mathbf{r}_{ij}}{x_m}\label{eqn:gradient}
\end{align}
%
The last two terms of the sum will be the same independent of the choice for the
potential function $\varepsilon(r_{ij})$. The implementation therefore 
focused on reducing redundancies by making use of class inheritance features. As
explained in the following paragraphs, this makes exchanging the type of
two-body potential trivial.

For this, a base class called \verb|pairPotential| was defined. It is an
abstract class and therefore cannot be instantiated. Its private fields hold
declarations of virtual methods for calculating energy $E(r)$, first derivative
$\dv*{E}{r}$ and second derivative $\dv*[2]{E}{r}$ based on particle distance
$r$. They are declared virtual, because they will be overwritten with the
respective functions in the derived classes of the actual potentials. In the
bare \verb|pairPotential| class those functions are only declared but never
defined and cannot be used for calculations. The public members of the class are
the constructor and the user-accessible functions for calculating energy,
gradient vector and Hessian matrix as well as the optimiser.

As an example, the \Cpp implementation for the member function that calculates
the energy of the system is shown in listing~\ref{code:energy}.
%
\begin{lstlisting}[caption={Implementation of the redundant part of the energy calculation.},label=code:energy,float=htb]
double pairPotential::calcEnergy (structure &S) {
  double f(0);
  for (int i = 0; i < S.nAtoms(); i++) {
    for (int j = i + 1; j < S.nAtoms(); j++) {
      f += this->E (coord3d::dist (S[i],S[j]));
    }
  }
return f;
}
\end{lstlisting}
%
The method takes a \verb|structure| object as input and uses the virtual energy
function to calculate the energy contributions of all unique pairs of spheres.
The method \texttt{nAtoms()} returns the number of spheres in the
\texttt{structure} and \texttt{dist} calculates the Euclidean distance between
two spheres \texttt{i} and \texttt{j}.

\begin{lstlisting}[caption={Implementation of the distance dependant energy for the Lennard-Jones potential.},label=code:LJenergy,float=htb]
class LJ : public pairPotential {
  private:
    double E (double distance) {
      return (_epsilon / (_exp1/_exp2 - 1)) 
             * ( (pow (_rm / distance, _exp1)) - (_exp1/_exp2) 
             * (pow (_rm / distance, _exp2)) );
    }
    //methods for first and second derivatives go here
};
\end{lstlisting}

To define a potential, a derived class that overwrites the virtual function
declarations is required. The virtual member functions of the base class are
overwritten in the private field of the derived class by providing properly
defined methods. Additionally, any parameters, that the potential form depends
on, are declared in the private fields. In the public fields, constructors for
the respective potential as well as a function that reads the potential
parameters from a user provided file need to be declared. The last function is
important as it also creates a pointer to the potential object on heap memory,
which is necessary, because the exact nature of the potential is not known at
compile time. For the \ac{LJ} potential the energy function can be defined as
shown in listing~\ref{code:LJenergy}. The variables starting with an underscore
are data members defined on object creation and refer to the two exponents
\texttt{\_exp1} and \texttt{\_exp2} the equilibrium distance \texttt{\_rm} and
the depth of the potential energy well \texttt{\_epsilon}. A different potential
can be implemented in the same way. A new derived class has to be defined,
containing the methods to calculate energy, first and second derivative.
Additionally, a function that reads parameters that define the respective
potential form has to be provided. By design, the function needs to return a
pointer to an instance of the class on heap memory. In case of the \ac{LJ}
potential this can be achieved as shown in listing~\ref{code:readPotential}.
%
\begin{lstlisting}[caption={Minimal example for the method \texttt{readPotential()}.},label=code:readPotential,float=htb]
LJ *LJ::readPotential () {
  //instructions to read parameters from file go here

  LJ *potential = new LJ(epsilon, rm, exp1, exp2);
  return potential;
}
\end{lstlisting}
%
Because memory on heap has to be managed by the user, the pointer should be used
in conjunction with \verb|unique_ptr| to ensure destruction of the object upon
exiting the scope. An example for this is given below.
%
\begin{lstlisting}
//read potentialName from settings file

std::unique_ptr< pairPotential > potential;
if (potentialName = "LJ") {
  potential.reset( LJ::readPotential() );
}
\end{lstlisting}

%\chapter{Computational Details}
%\label{sec:ComputationalDetails}
%
%The following sections detail the methodologies used in part~\ref{sec:results}.
%They have been altered slightly from the original
%publications\autocite{Trombach_HollowGoldCages_2016,%Trombach_stickyhardsphereLennardJonestypeclusters_2018,%Trombach_GregoryNewtonproblemkissing_2018}
%to match the style of this thesis.
%
%\section[Computational Details for
%Chapter~\ref{sec:goldendualfullerenes}]{Computational Details for
%Chapter~\ref{sec:goldendualfullerenes}\footnote{This section has been published
%previously in
%\citetitle*{Trombach_HollowGoldCages_2016}\autocite{Trombach_HollowGoldCages_2016}
%and is reproduced with permission from the publisher \textcopyright 2016
%WILEY-VCH Verlag GmbH \& Co. KGaA.}}
%\label{sec:ComputationalDetailsforGoldenDualFullerenes}



%\section[Computational Details for
%Chapter~\ref{sec:fromstickyhardspheretoLJtypeclusters}]{Computational Details
%for Chapter~\ref{sec:fromstickyhardspheretoLJtypeclusters}\footnote{This
%section has been published previously in
%\citetitle*{Trombach_stickyhardsphereLennardJonestypeclusters_2018}\autocite%{Trombach_stickyhardsphereLennardJonestypeclusters_2018}
%and is reproduced with permission from the publisher \textcopyright 2018
%American Physical Society.} }



%\section[Computational Details for
%Chapter~\ref{sec:thegregorynewtonclusters}]{Computational Details for
%Chapter~\ref{sec:thegregorynewtonclusters}\footnote{This section has been
%published previously in
%\citetitle*{Trombach_GregoryNewtonproblemkissing_2018}\autocite%{Trombach_GregoryNewtonproblemkissing_2018}
%and is reproduced with permission from the publisher \textcopyright 2018
%American Physical Society.}}