Search

The accurate prediction of protein-ligand binding free energies with tractable computational methods has the potential to revolutionize drug discovery. Modelling the protein-ligand interaction at a quantum mechanical level, instead of relying on empirical classical mechanical methods, is an essential step toward this goal. In this body of research, we explore the QM-PBSA method to calculate quantum mechanical free energies of binding based on full-protein linear-scaling density functional theory calculations. We apply the QM-PBSA method to the T4-lysozyme L99A/M102Q protein and investigate the convergence, precision, and reproducibility of the QM-PBSA method. Additionally, we compare three different exchange-correlation functionals and different empirical dispersion corrections. Building on our findings in the well-characterized T4-lysozyme we calculate quantum mechanical protein-ligand free energies of binding for a set of ligands binding to the pharmaceutically highly relevant bromodomain containing protein 4 (BRD4) after an extensive investigation of the protein system at the classical mechanical level. BRD4 plays a key role in many cancers. The inhibition of BRD4 can suppress the cancer growth of acute myeloid leukemia, diffuse large B cell lymphoma, prostate cancer, and breast cancer. We demonstrate the predictive power of QM-PBSA in BRD4 as compared to its classical mechanical analog MM-PBSA and show, beyond doubt, that full-protein quantum mechanical calculations are both viable and tractable on modern supercomputers and in an academic context.

Nanopore DNA sequencing is a well-established technology that has accelerated advancements in many fields, including medical research. Over the years, research has focussed on optimising protein nanopores for DNA sequencing. Optimisation strategies broadly focus on (1) slowing the translocation of DNA to increase the time available for base recognition and (2) improving the resolution of detection to attain single-base sensing. Molecular dynamics (MD) simulations have been invaluable in obtaining molecular-level insights to pave the way for informed nanopore optimisation. In this thesis, MD simulations were used to study DNA translocation through nanopores and elucidate the design principles for optimising nanopores for DNA sequencing. In the first chapter, the translocation of short and longer single-stranded (ss)DNAs was studied through protein-inspired hydrophobic nanopores with dual-constrictions. It was found that DNA translocation is slowed down by aromatic residues, and when combined with a narrow geometry, DNA retains a largely linear conformation during translocation and without forming secondary structures that can impede DNA sequencing. Following this, the proteins CsgG and the CsgG-CsgF complex were characterised in terms of their conformational dynamics and ability to allow DNA translocation. Eyelet loops forming the CsgG constriction were found to exhibit large variations in their mobility, with at least one loop moving upwards into the vestibule under an applied electric field. CsgF was found to stabilise CsgG and the eyelet loop region. Subsequently, the translocation of short ssDNA through CsgG and the CsgG-CsgF complex was studied. The speed of DNA translocation was found to be primarily influenced by DNA interacting with key residues in the CsgG constriction region. DNA is retained in a more linear conformation during translocation through the dual-constriction hydrophobic channel formed by the CsgG-CsgF complex compared to CsgG. Next, the translocation of longer ssDNA in an applied electric field was studied through these proteins/protein complex. These simulations revealed that the eyelet loops of the CsgG constriction region are mobile during DNA translocation, and the stochastic nature of their mobility perturbs the pore geometry which may give rise to noise in the ionic current through the nanopores. In the last chapter, Markov State Model methodology was employed to characterise the kinetics of the mobility of the CsgG eyelet loops under an applied electric field. The model construction was limited by the duration of the MD simulations. The data and analyses presented in this, and previous, chapters emphasise the need for a model that describes the complex conformational dynamics of the CsgG eyelet loops.

The behaviour of lubricants at operational conditions, such as at high temperatures and pressures, is a topic of great industrial interest. In particular, viscosity and the viscosity-pressure relation are especially important for applications and their determination by computational simulations is very desirable. In this thesis we evaluate methods to compute these quantities based on fully atomistic molecular dynamics simulations which are computationally demanding but also have the potential to be most accurate. We tested several molecules that are used as lubricants, such as 9,10-dimethyloctadecane, main component of PAO-2 base oil, which was used as the main lubricant for our simulations. The methods used for the viscosity simulations are the Green-Kubo equilibrium molecular dynamics (EMD-GK), the direct computation of viscosity from shear during non-equilibrium MD (NEMD) and the use of confined NEMD, where the fluid is confined within explicitly defined iron oxide wall surfaces, at pressures of up to 1.0 GPa and various temperatures (40-150 degrees Celsius). We present the theory behind these methods and investigate how the simulation parameters affect the results obtained, to ensure viscosity convergence with respect to the simulation intervals and all other parameters. We show that by using each method in its regime of applicability, we can achieve good agreement between simulated and measured values. NEMD simulations at high pressures captured zero shear viscosity successfully, while at 40 degrees Celsius EMD-GK is only applicable to pressures up to 0.3 GPa, where the viscosity is lower. In NEMD, longer and multiply repeated simulations reduce the standard deviation of viscosity, which is essential at lower pressures. Additionally, by using confined NEMD simulations, it was demonstrated that the film thickness of the fluid affects viscosity, and as we increase the number of lubricant molecules, we approach the viscosity value of the bulk fluid derived from NEMD simulations. Another aspect of these methods is the choice of the utilised force field for the atomic interactions. This was investigated by selecting three different commonly used force fields. We have explored several methods for calculating viscosity and we obtained results of particular industrial interest.

Molecular systems with functional domains serve as a practical motivation for understanding the factors that contribute to the interaction energy. The ability to decompose the interaction energy of a group of interacting subsystems is an important method in studying the chemical nature of the interactions. Energy decomposition analysis (EDA) is a family of schemes that allows such dissection of the interaction energy into chemically relevant components depending on the scheme used. Since different EDA schemes decompose the interaction energy differently, the interpretation of the resulting components differs among schemes. However, various EDA schemes provide complementary insights into the interactions between chemical entities. In this work, two EDA schemes are developed or extended: Hybrid Absolutely Localized Molecular Orbitals (HALMO) and Combined Localized Molecular Orbitals (CLMO). Both EDA schemes have been implemented as part of a linkable library alongside the computational chemistry package, ONETEP. Since ONETEP is a linear-scaling software package, an important application of such decomposition analysis is in the study of large, nontrivial molecular systems for more insightful understanding of chemical interactions, which in turn can lead to more accurate and focused design of chemical systems. Systems such as biomolecules usually contain several self-stabilizing domains that can fold independently and have important functions. Defining the fragments of a supermolecule is necessary in EDA, and if done appropriately given the context of a particular application, the fragmentation of a biomolecule can elucidate the intramolecular interactions that contribute to the functions of the system as a whole. Two major methods of self-consistent field for molecular interaction (SCF MI) are examined and made more mathematically transparent. SCF MI was originally designed to exclude basis set superposition error (BSSE) from molecular interactions. However, SCF MI is also used in separating charge transfer from polarization in HALMO EDA. Studying several SCF MI methods provides HALMO EDA alternatives for separating charge transfer as an EDA component. Due to ONETEP's linear scaling, large biomolecules or large samples of biomolecules could be studied for their interactions or distributions of specific properties. The primary type of biomolecule studied in this work is double-stranded DNA (dsDNA) and how its stability relates to guanine-cytosine (GC) content, which is a measure of the amount of guanines or cytosines in nucleic acids. By examining the interaction energies in terms of EDA components, contributions to the variabilities within and across GC-content groups are examined and are correlated with differing stabilities despite having same GC content. Ensemble density-functional theory (EDFT) optimizes the molecular orbitals of a system at finite electronic temperatures and allows occupancies to be fractional when the band gap is sufficiently small. EDA methods are normally developed for the pure state and pose difficulties in decomposing the interaction energy of a conductor due to the fractional occupancies that are part of the optimization process in EDFT. As such, fractional occupancies have been incorporated in the EDA optimization process of the fragmented species, thereby allowing EDA to be applied to systems of relevance to catalysis and metallic systems. The adaptations of EDA and SCF MI to metallic systems are novel and were validated using samples from catalysis and batteries, and HALMO EDA has provided reasonable decompositions of interaction energies and revealed some trends from SCF MI that correlate with charge distributions and chemical intuition.

Understanding solvent effects at the metallic/liquid interface is critical to improving and analysing heterogeneous catalytic processes. In addition to a growing body of experimental work, computational studies are elucidating how the presence of water affects both the electronic structure and the adsorption thermodynamics of the metallic surface. However, computational methods such as ab initio Molecular Dynamics (AIMD) require extensive configurational sampling to obtain equilibrated thermodynamic quantities, precluding their use for wide-ranging studies. In contrast, by considering the dynamic degrees of freedom as an average, implicit solvent methods provide a route to tractable computational simulations of the aqueous environment, while maintaining a quantum mechanical description of metallic/adsorbate interactions. This thesis, in collaboration with the Pacific Northwest National Laboratory (PNNL), describes how implicit solvent approaches can be applied as an inexpensive method of evaluating both the electronic structure of the metal/liquid interface, and the free energy change of adsorption in the aqueous phase for a range of organic adsorbates. To ensure these calculations are performed accurately and efficiently, developments were made to the linear scaling Density Functional Theory (DFT) code, ONETEP as part of this work. These developments include the implementation of the soft sphere dielectric cavity model, which gives the flexibility to parameterize the solvent model for individual atomic centres. This contrasts with the original electron density based cavity model, which applies a global cavity parameter, leading to poor descriptions of the free energy changes of solvation for systems with mixtures of light organic and heavy metallic species. A surface accessible volume term was also implemented for the non-polar solvation term, which improves the correlation with experimental solvation free energies compared to the surface area non-polar term. Furthermore, a Pulay Hamiltonian mixing routine was implemented in the Ensemble DFT (EDFT) scheme of ONETEP. This approach confers significantly improved convergence behaviour for single point energy calculations performed in this work. This enables more efficient and accurate simulations of metallic systems, allowing for the evaluation of larger systems studied in later chapters. Utilising the implemented soft sphere model, this work assesses the ability of the soft sphere model to capture the potential of zero charge and the work function of the metallic/liquid interface. By reparameterizing the implicit solvent model in terms of the work function values calculated from snapshots of an AIMD simulation, we were able to capture the salient electronic structure changes of the solvated metallic surface and electrochemical properties. Then, combining the accelerated EDFT scheme and the implicit solvent parameterization method used for the potential of zero charge, this thesis concludes with a continuum solvent approach for calculating the aqueous phase adsorption free energy of organic molecules to the Pt(111) surface. In this work, approximations are derived for the entropies of solvation for the metallic surface based on analytical statistical thermodynamic expressions. These approximations allow us to parameterize the implicit solvent mode ∆G_solv for the metallic surface, enabling adsorption free energy with reasonable accuracy for a range of coverages and orientations. This opens a route for computationally inexpensive evaluations of adsorption processes at the aqueous Pt(111) interface, which can provide an atomistic understanding of adsorption processes in support of experimental studies. The work presented in this thesis shows the usefulness of the implicit solvent method in studies of heterogeneous catalytic processes and electrochemical interfaces. The techniques described in this work show that thermodynamic and electrochemical properties can be calculated in a computationally tractable manner with implicit solvent. In future, this could enable high throughput studies for a range of metallic surfaces and adsorbates, aiding the design of catalysts for a range of applications.

Heterogeneous catalysts are used in a huge number of applications across many industries. Determining how each of a multitude of factors influence catalytic activity is a complex and difficult task, but is necessary if new, improved catalysts are to be developed. While much of this work is applicable to heterogeneous catalysis in general, the oxygen reduction reaction (ORR) that takes place in hydrogen fuel cells has been a central theme. This thesis describes how computational studies can be used to build on our understanding of such catalysts, focusing on the effects of nanoparticle morphology, alloying and those of the support. In order to capture these effects, ever more realistic simulations must be performed. We compare real platinum nanoparticle structures from microscopy experiments with simulated models, and find that the richness of morphologies seen in the real nanoparticles is important for activity. A method to build comparable structures via annealing is also described. For cobalt platinum alloyed nanoparticles, which are relevant to the ORR, a direct comparison of structures is much more challenging. We have compared measurements of structural strain between theory and experiment in order to help determine the cobalt distribution within such nanoparticles, which remains an important and unresolved question, and is particularly important for alloyed nanoparticle stability. Finally, the effects of both titania and alumina supports are shown to be significant for ligand adsorption, and models should include them where possible. Titania is shown to reduce the binding strength of oxygen and carbon monoxide when used as a support for platinum nanoparticles. These studies represent different steps in the direction of realistic simulations, where more factors are accounted for than in current routine calculations. The goal of any simulation is to capture all the important details of a real system in order to better reproduce the results, and ultimately make consistent and accurate predictions. Such simulations would be invaluable to catalyst design, reducing the reliance on expensive trials of real systems.

Organic photovoltaic (OPV) devices rely on the mixture between a conjugated copolymer (electron donor) and an electron acceptor material (typically, but not necessarily, functionalised fullerenes): this active layer is known as the bulk heterojunction, and it is crucial for the device operation, as this is where excitons are split into free electrons and holes to produce current. A deep understanding of the role of the molecular structure of materials on the device physics is necessary to achieve better performances, and to this end, computer simulations are undoubtedly a powerful tool. However, the bulk heterojunction is a complex system, and for theoretical models to be representative, these should be composed of thousands of atoms, a size out of the reach of atomistic quantum mechanics simulations. To overcome this limitation, this project involved the use of the ONETEP code, which, thanks to its linear-scaling computational cost with respect to the system size, allowed to carry out ab initio calculations, within the density functional theory (DFT) framework, on OPV materials and models of bulk heterojunctions on a far larger scale than possible before. Regardless of the device morphology and architecture, fundamental for the exciton splitting is the energy level alignment of the donor and the acceptor components, and several ways exist to fine-tune the electronic properties of these materials. Nevertheless, is it possible to find novel and alternative routes to the well-known strategies currently employed in the laboratory? As for the donor polymer, the results here presented suggest that acting on the polymerisation statistics, that is, the ratio of different blocks in the polymer chain, significantly affects the electronic structure of such materials, with changes in the band gap of the same order of magnitude of those induced by the widely used functionalisation approach. The acceptor fullerenes, on the iii other hand, are generally more challenging to functionalise, and consequently their electronic structure cannot be trivially tuned. However, one could ideally circumvent the issue via the intercalation of different solvent molecules in the crystal phase of fullerenes, in order to attempt to indirectly modify the energy of the frontier orbitals and the band gap. Indeed, results highlighted the crucial role of the solvent in modifying both the electronic and the optical properties of solid-state fullerenes through the formation of fullerene-solvent ⇡-⇡ interactions, which disrupt the close packing of solvent-free fullerenes. Interestingly, more appreciable changes were observed in the properties of pure rather than functionalised fullerenes. Another, yet fundamental, aspect of OPV explored here is the polymeracceptor interaction in the bulk heterojunction, with the acceptor material consisting of fullerene and the more recently introduced non-fullerene acceptors (NFAs). Although attempts to model the polymer-fullerene interface to investigate its excited-state properties are numerous, these have been limited to, within the framework of ab initio atomistic simulations, small 1-to-1 short-oligomer-fullerene pair models. For the first time, DFT calculations for ground and excited state were performed on model interfaces of realistic size composed of more than a single polymer chain and dozens of fullerene molecules, allowing to gain new insights into the physics of exciton generation and splitting. For instance, it was observed that the probability of charge transfer to occur is deeply influenced by the polymer block statistics, and that exciton dissociation is favoured by large polymer phases rather than large fullerene phases, although these are still beneficial. Evidence of long-range charge-transfer states in the low-energy part of the excited-state spectrum was also observed. On the other hand, models of polymer-NFA interfaces are still scarce in the literature, as NFA phases are more intrinsically complex to model than fullerene phases. Critical for both charge mobility and device performance is the NFA solid-state arrangement with respect to the polymer. By constructing large polymer-NFA model interfaces it was possible to highlight and confirm the importance of intermolecular ⇡-⇡ stacking interactions in the NFA phase, as it was found that these deeply influence the exciton delocalisation, the exciton splitting rate, and the mobility anisotropy. This work, which is the outcome of collaboration with Merck, was enabled by the linear-scaling capabilities of the ONETEP code, which allowed to study large-scale realistic models of OPV. This thesis provided novel and important insights into different aspects of organic photovoltaics, both in terms of material design and device physics.

Understanding the factors controlling a catalyst activity, selectivity, and stability is a complicated task with potential applications for a large number of technologies with high impact for our society. For example, optimised catalysts can have a central role to improve the durability, efficiency and decrease the cost of technologies such as fuel cells, which are devices able to convert the chemical energy from molecules into electricity by performing chemical reactions that are specific to the used fuel and the type of fuel cell. Fuel cell catalysts are commonly made from metallic nanoparticles on different types of supports, meaning that characteristics such as the nanoparticle composition, size, shape, composition of the support, electrolyte and many others can be simultaneously controlled to tune the catalyst towards a specific goal. The complexity involved in the catalyst optimisation makes the problem extremely exciting and challenging, requiring significant effort from many different research areas ranging from synthetic chemistry and electrochemistry to computational chemistry. This thesis describes computational studies on metallic nanoparticles, focusing on nanoparticle size effects and its interplay with other variables that are important in the context of fuel cell catalysts such as the presence of support and adsorbate coverage effects. We also present a framework in ONETEP's linear-scaling DFT formalism for the implementation of local and angular momentum projected density of states (l-p-DOS), which is an important tool to study metallic nanoparticles used as catalysts. The four different bases used to project the density of states are tested, and its results are compared against other DFT code, helping to validate our l-p-DOS implementation. The results obtained for metallic nanoparticles show similar trends with all the implemented options, demonstrating the reliability of the method for such studies. The ONETEP code with the implemented l-p-DOS functionality is used to perform a set of large-scale DFT calculations to study Pt nanoparticles isolated and supported on pristine graphene. The results show a weak metal-support interaction, with the adhesion energy per Pt atom decreasing with the nanoparticle size and being dominated by dispersion interactions for larger nanoparticles. The interaction with the support induces geometric and electronic changes in the nanoparticle, with the inter-atomic distances expanding (contracting) for Pt facets close (far) to the support and with charge redistribution happening at the interface between Pt clusters and graphene. The changes in the geometric and electronic properties induced by the interaction with the support are size dependent, correlate with the interaction strength between cluster and support, and generate shifts in the d-band centres of the Pt nanoparticles. The isolated and supported nanoparticles were used to study how the nanoparticle size and presence of support alters the interaction of the nanoparticles with atomic oxygen, carbon monoxide, and ethanol. We show that nanoparticle size decreases strengthen the adsorption energies of O and CO, which can hinder the activity of these nanoparticles towards important reactions in the context of fuel cells. The size effect is observed to control the adsorption energies for O and CO to a larger extent than for ethanol. We also note that the presence of graphene weakens the adsorption energies for all adsorbates, with this effect being more significant for smaller nanoparticles. Finally, we study how the adsorbate coverage changes the adsorption of atomic oxygen on Pt nanoparticles and how coverage effects vary with the nanoparticle size. We observe that the increase in O coverage weakens the adsorption energies per O atom, with this effect being more significant for the larger simulated nanoparticles. These studies can help to understand different influences that the nanoparticle size can have by simulating non-trivial combinations of size and support effects and size and coverage effects which could have implications to designing more efficient catalysts for fuel cells and other applications.

This thesis will focus on DFT for calculations of large metallic nanoparticles. It will show new algorithms that were developed for reduced scaling DFT methods for metals; the testing, verification and design of new descriptors for predicting the catalytic activity of metallic nanoparticles; application of large-scale DFT calculations to model nanoparticle sequences to show size and oxygen adsorption coverage trends, and finally the application of these techniques and knowledge to perform a study of oxygen adsorption on real-world, experimentally determined platinum nanoparticles in collaboration with the Nellist group at Oxford materials. We explore the binding of atomic oxygen to cuboctahedral platinum nanoparticles of up to 1000 atoms using DFT calculations in ONETEP. We demonstrate convergence to the infinite slab limit for single oxygen adsorption in chapter 4 and correlate adsorption strength against popular descriptors for catalytic activity, such as the d-band centre approach. This approach is possible because of work which will be described in chapter 3 to implement angular momentum projected density of states calculations in ONETEP. The effects of oxygen coverage on the Pt55 and Pt147 cuboctahedral nanoparticles will also be analysed, which serves to advance our simulations towards realistic conditions. We show in our investigation into half monolayer, hemispherical oxygen coverage on platinum nanoparticles that oxygen tends to gravitate towards the edges and lower coordinated sites in the nanoparticle and away from the centres of facets. This effect correlates with the site specific, single oxygen adsorption energies on Pt309 and experimental platinum nanoparticles which is presented in chapter 5. We show that when subdividing the binding of monolayers of oxygen into only (111) and (100) facets that these have a lower adsorption strength per oxygen atom than combined (100) and (111) facets as well as lower binding strength than single oxygen adsoprtion. In the next part of the study, which is discussed in chapter 5, we show large scale DFT calculations on real platinum nanoparticles, which were measured by the Nellist group at Oxford materials using advanced electron microscopy techniques. These DFT calculations provide the electronic structure of the experimentally measured nanoparticles, which allowed us to apply electron density based catalytic activity descriptors to the nanoparticles, such as the d-band centre approach, or our own electronic density based descriptor described in chapter 3. We find that surface roughness of the experimental nanoparticles contributes to more potential oxygen binding sites with low electron density, which correlatates with stronger oxygen adsorption strength in our model, when compared with the relative smoothness of cuboctahedral and truncated octahedral facets. In the analysis which is presented in chapter 5, the proportion of sites which lie within 0.2 eV of the oxygen binding strength required for optimum catalytic activity is predicted with high efficiency, based on our catalytic activity descriptor. Finally, in chapter 6 we describe a new method for large scale DFT calculations on metallic systems which we call the AQuA-FOE method. We show how this method can have a computational cost which increases effectively linearly with the number of atoms. The AQuA-FOE method works by implicitly heating and quenching the electrons in the system to find the oneparticle density matrix, while conserving the electron number. We show validation of this method inside the EDFT procedure by comparing numerically with the diagonalisation based EDFT that is already implemented in ONETEP showing agreement in the energies to better than 10⁻⁵ EH per atom. We will also demonstrate the effectively linear-scaling computational cost of our method with calculation times on regular truncated octahedral Palladium nanoparticles ranging from 2,406 to 12,934 atoms.

Nowadays, Kohm-Sham density functional theory (KS-DFT) calculations are routinely employed in several research fields, due to the ability of KS-DFT to provide great accuracy for a wide class of molecular systems and materials. Unfortunately, conventional KS-DFT calculations, although very powerful, require a computational cost that goes with the cube of the system size, also known as ON³ scaling, undermining in practice the study of large systems. The advent of linear-scaling, or O (N ), DFT (LS-DFT) methods, relying on the locality of the electronic matter, has enabled calculations on increasingly large systems, up to tens of thousands of atoms. The central tenet of linear-scaling methods is the exponential decay in real space of the single-particle reduced density matrix. This property allows to enforce localisation constraints on the electronic structure, significantly reducing the size of the matrices, such as the Hamiltonian matrix, and increasing their sparsity. The single-particle density matrix in the LS-DFT formalism is expanded in terms of a set of atom-centred, strictly localised functions. Employing periodic boundary conditions (PBCs), the energy is minimised with respect to all the degrees of freedom in the density matrix, which allows to attain chemical accuracy using a high-resolution minimal basis set. The combination of localisation constraints and sparse algebra form the substrate for O (N ) calculations. In this thesis, we used Onetep, a linear-scaling DFT program, to carry out our calculations. The aim of our research is to combine molecular dynamics simulations, within the Born-Oppenheimer approximation (BOMD), with linear-scaling DFT methods. In particular, our main goal is to advance current methodology by developing new algorithms to better exploit locality in BOMD and to reduce the computational load while maintaining DFT accuracy. Dipole moment autocorrelation functions can directly be employed to obtain the IR spectrum of a given molecular system. DFT-MD simulations offer the perfect tool to generate accurate autocorrelation functions which automatically take into account the anharmonicity of the potential energy surface and temperature effects. Computational IR spectroscopy plays a pivotal role in the understanding of conformational changes of biomolecules, which tend to show several almost-degenerate conformers at room temperatures (floppy molecules). It is particularly valuable when interpreting their fingerprint in solution in combination with experimental spectra. We have implemented two algorithms for the computation of the local electronic dipole moments of molecules in solution. Both methods demand a strategy to partition the density. These methods enable the computation of IR spectra of large molecules, such as polypeptide, in explicit solvent. In the resulting IR spectrum the effect of the solvent on the target molecule is automatically captured, whereas its IR signature is removed. We expect these new functionalities to be very helpful in the understanding of how bio-molecules interact with the solvent at room temperature and the effect of these interactions on conformational changes. Computationally, the most demanding step in molecular dynamics simulations is the evaluation of energies and forces. This has particular severe consequences on BOMD-based approaches. In fact, a self-consistent field (SCF) step is required at each MD step, which in turn requires multiple energy evaluations. As a consequence, the SCF loop has a major effect on the computational load and overall wall-time. MD Schemes that are capable to by-pass the SCF loop altogether, e.g. Car-Parrinello MD or fixed charge force fields in classical MD, are inherently faster in terms of wall-time per MD step, although they usually demand a much smaller time-step. Moreover, the quality of the converged density matrix is crucial for energy conservation and forces in the LS-DFT BOMD approaches. In theory, the self-consistent solution does not depend on the initial guess. In practice, the SCF optimisation is always incomplete, leading to memory effects and the breaking of time-reversal symmetry, which gives rise to systematic errors in energy gradients that manifest as a drift in microcanonical energy. To ameliorate this problem, we present two integration schemes based on an extended-Lagrangian (XL) approach which introduces extended or auxiliary electronic degrees of freedom to generate good quality time-reversible initial guesses in the SCF loop. Both schemes are improvements over the original XL formulation, which suffered from numerical noise accumulation. The first approach, known as dissipative XL (dXL), introduces a dissipative-like term in the Verlet integration (modified Verlet) of the auxiliary degrees of freedom; the second approach, known as inertial XL (iXL), introduces a thermostat, hence requiring a velocity Verlet integrator. We have implemented both schemes in Onetep and studied their performance using liquid water as test case. In collaboration with the authors of the schemes, we have analysed the energy driftt in both classical polarisable force field MD calculations and LS-DFT BOMD. We have found that both schemes are very efficient in reducing the number of SCF iterations while maintaining good energy driftss and have similar performance. We believe that our implementation and analysis will be very beneficial for future applications both with LS-DFT BOMD and classical polarisable force field MD since both schemes constitute important algorithmic improvements that markedly extend the timescales accessible to classical and LS-DFT MD simulations alike.

Kohn-Sham density functional theory (DFT) is an extraordinarily powerful and versatile tool for calculating the properties of materials. In its conventional form, this approach scales cubically with the size of the system under study. This scaling becomes prohibitive when investigating larger arrangements such as biomolecules and nanostructures. More recently linear-scaling approaches have been developed that overcome this limitation, allowing calculations to be performed on systems many thousands of atoms in size. An example of such an approach is the ONETEP code which uses a plane wave-like basis set and is based upon the use of spherically-localised orbitals. A simple yet common calculation performed using ab initio codes is the total (ground state) energy calculation. By comparing the energy of isolated parts of a system to the energy of the combined system, we are able to obtain the energy of interaction. This quantity is useful as it provides a relative measure of the enthalpic stability of an interaction which can be compared to other systems. Equally, however, this quantity gives little indication of the driving forces that lead to the interaction energy we observe. A number of approaches have been developed that aim to identify these driving forces. Energy decomposition analysis (EDA) refers to the set of methods that decompose the interaction energy into physically relevant energy components which add to the full interaction energy. Few studies have applied EDA approaches to larger systems in the thousand-atom regime, with the vast majority of investigations focussing on small system studies (less than 100 atoms in size). These methods have shown varying degrees of success. In this work, we have evaluated the suitability of a selection of popular EDA methods in decomposing the interaction energies of small biomolecule-like systems. Based on the results of this review, we developed a linear-scaling EDA approach in the ONETEP code that separates the intermolecular interaction energy into chemically distinct components (electrostatic, exchange, correlation, Pauli repulsion, polarisation, and charge transfer). The intermediate state used to calculate polarisation, also known as the absolutely localised molecular orbital (ALMO) state, has the key advantage of being fully antisymmetric and variationally optimised. The linear-scaling capability of the scheme is based on use of an adaptive purification approach and sparse matrix equations. We demonstrate the accuracy of this approach in reproducing the energy component values of its Gaussian basis counterparts, and present a remedy to the limitation of polarisation and charge transfer basis set dependence that is based on the property of strict localisation of the ONETEP orbitals. Additionally, we show the method to have mild exchange-correlation functional and atomic coordinate dependence. We have demonstrated the high value of our method by applying it to the thrombin protein interacting with a number of small binders. Here, we used our scheme in combination with electron density difference (EDD) plots to identify the key protein and ligand regions that contribute to polarisation and charge transfer. In our studies, we assessed convergence of the EDA components with protein truncation up to a total system size of 4975 atoms. Additionally, we applied our EDA to binders that had been partitioned into smaller fragments. Here, we accurately quantified the bonding contributions of key ligand moieties with particular regions of the protein cavity. We assessed how accurately the ligand binding components are reproduced by the fragment contributions using an additivity measure. Using this measure, we showed the fragment binding components to add up to the full ligand binding component with overall minimal additivity error. We also investigated the energy components of a series of small thrombin S1 pocket binders all less than 30 atoms in size. In this study, we demonstrate the EDA and EDD plots as tools for understanding the relative importance of different binder structural features and positionings within the pocket. Overall, we show our EDA method to be a stable and powerful approach for the analysis of interaction energies in systems of large size. The application of this method is not limited to biomolecular studies, and we expect that this approach can be readily applied to analyses within other fields, for example materials, catalysts, and nanostructures.

The work presented within this thesis uses quantum mechanical (QM) calculations to improve free energies of binding computed with classical (MM) force fields. Initially a direct approach was taken, where snapshots were taken at equally spaced distances throughout the classical simulation and each structure underwent a quantum single point energy calculation. This direct approach was possible by using the Zwanzig equation. However, one disadvantage of using the Zwanzig equation is it's extreme sensitivity to fluctuations in the energy difference. This led to the quantum corrected free energy being dominated by a few snapshots and convergence could not be achieved. This led to the application of an acceptance criterion, where instead of just using each evenly spaced structure from the classical simulation, each structure would have to be accepted into a target potential. In the work presented here, our target potential was a QM potential. For a simple test system of N2 in vacuum we achieved a high acceptance and converged free energies, however, for more complex systems little to no acceptance was found. The poor acceptance can be attributed to the difference between the MM and QM potential energy surfaces. Similarities were found, however, on the minima of these potential energy surfaces between the MM and QM, which led to the application of a bias to ensure that sampling was only taken from the minima. However, similar to the Zwanzig equation, this method proved to be too sensitive to difference in energy, thus convergence could not be achieved. In order to "smooth" the transition between the MM and QM a "stepping stone" approach was used. The first step was to accept structures from a classical simulation to a QM/MM ensemble, then we used a direct approach again using the Zwanzig equation to move from the QM/MM potential to the QM. Using this approach, we find very small convergence errors (< 1 kJ/mol). This method was validated by calculating hydration free energies for a variety of ligands. Finally, the free energy of binding was calculated for trypsin with several benzamidine derivatives using a QM-PBSA approach, which involved running QM calculations on the entire protein-ligand complex. The final results, however, showed no overall improvement of the calculated free energies between the MM and QM. It was found that the inclusion of QM methods lowered the free energy in each case.

In the drug discovery process, there are multiple factors that make a successful candidate other than whether it antagonises a chosen active site, or performs allosteric regulation. Each test candidate is profiled by its absorption into the bloodstream, distribution throughout the organism, its products of metabolism, method of excretion, and overall toxicity; summarised as ADMET. There are currently methods to calculate and predict such properties, but the majority of these involve rule-based, empirical approaches that run the risk of lacking accuracy as one's search of chemical space ventures into the more novel. The lack of experimental data on organometallic systems also means that some of these methods refuse to predict properties on them outright, losing the opportunity to exploit this relatively untapped area that holds promise for new antibacterial and antineoplastic pharmaceutical compounds. Using the more transferable and definitive quantum mechanical (QM) approach to drug discovery is desirable, but the computational cost of conventional Hartree-Fock (HF) and Density Functional Theory (DFT) calculations are too high. Using the linear-scaling DFT program, onetep, we aim to exploit the benefits of DFT in calculations with much larger fragments of, and in some cases entire biomolecules, in order to demonstrate calculations which could ultimately be used in developing more accurate methods of profiling drug candidates, with a computational cost that albeit still high, is now feasible with the provision of modern supercomputers. In this thesis, we first use linear-scaling DFT methods to address the lack of electron polarisation and charge transfer effects in energy calculations using a molecular mechanics forcefield. Multiple DFT calculations are performed on molecular dynamics(MD) snapshots of small molecules in a waterbox, with the aim of computing a MM!QM correction term, which can be applied to a forcefield binding free energy approach (such as thermodynamic integration) which will process a far greater number of MD snapshots. As a result, one will obtain the precision from processing very large numbers of MD snapshots of biomolecular systems, but the accuracy of QM. To improve efficiency of the QM phase of the overall method, we use electrostatic embedding to model the regions of the waterbox that are far from the solute, yet are still important to include. As this is a relatively new module in onetep, we present validation data prior to its use in the main work. Secondly, we validate different methods of calculating the pKa of a wide variety of molecules: from small, organic compounds, to the organometallic cisplatin, with the ultimate goal being of such calculations to eventually address questions such as, assuming oral intake, where in the gastrointestinal tract will a drug molecule be absorbed into the bloodstream, and how much of the original dose will be absorbed. These calculations are then scaled up significantly to examine the potential of using linear-scaling DFT to calculate the pKa of specific residues in proteins. This is performed with a 305-atom tryptophan cage, the 814-atom Ovomucoid Silver Pheasant Third Domain(OMSVP3) and a 2346-atom section of the T99A/M102Q T4-lysozyme mutant. We also highlight the challenges in calculating protein pKa. Finally, we study the hydrogen-abstraction reaction between cyclohexene and cytochrome P450cam, through onetep single point energy calculations of a 10-snapshot adiabaticreaction profile generated by the Mulholland Group(University of Bristol). Following this, the LST and QST methods of determining the transition state (available through onetep) are used, with the aims of determining the importance of the protein surrounding the active site in regards to the activation energy and structural geometry of the calculated transition state. The LST and QST methods are also validated, through modelling of the SN2 reaction between fluoride and chloromethane. The aim of this part of our work is to eventually assist in developing a metabolism (and toxicity) model of the different isoforms of cytochrome P450. Overall, this thesis aims to highlight not only the capability of linear-scaling DFT in becoming an important part of biomolecular simulation, but also the challenges that one will face upon scaling up calculations that were previously simple to perform, based on the small size of the system being modelled.

Kahn-Sham density functional theory (OFT) calculations yield reliable accuracy in a wide variety of molecules and materials. The advent of linear-scaling OFT methods, based on locality of the electronic matter, has enabled calculations on systems with tens of thousands of atoms. Localisation constraints are imposed by expanding the Kahn-Sham states in terms of a set of atom-centred, spherically-localised functions. Chemical accuracy is then achieved via a self-consistent optimisation using a high-resolution basis set. This formalism reduces the size of, and brings predictable sparsity patterns to, the matrices expressed in this representation, such as the Hamiltonian matrix. In this work, we used the ONETEP program for DFT calculations, which is based on the abovementioncd principles. The vision behind our research is to advance the method by developing new and robust algorithms to enable novel applications based on localised orbitals.

The accurate prediction of protein-drug binding affinities is a major aim of computational drug optimisation and development. A quantitative measure of binding affinity is provided by the free energy of binding, and such calculations typically require extensive configurational sampling of entities such as proteins with thousands of atoms. Current binding free energy methods use force fields to perform the configurational sampling and to compute interaction energies. Due to the empirical nature of force fields and the neglect of electrons, electron polarisation and charge transfer are not accounted for explicitly. This can limit the accuracy with which interactions are calculated and consequently the free energies obtained. Ideally ab initio quantum chemistry approaches should be used as these explicitly include the electrons. However, conventional ab initio approaches are not suitable due to their prohibitively high computational cost and unfavourable scaling. In this thesis we use large-scale ab initio quantum chemistry calculations within the Density Functional Theory (DFT) method to address the above mentioned limitations of force fields. To obtain quantitative results with ab initio approaches it is important to converge the calculations with the size of the basis set. For this reason we have used the ONETEP program, which is capable of linear-scaling DFT with near-complete basis set accuracy. A well known binding free energy approach is the Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA), which obtains free energies from evaluation of the energy of configurations in an implicit solvent model. We present the first application of a “QM-PBSA” approach to a protein-ligand system containing over 2600 atoms. In this QM-PBSA approach the energies of the configurations in vacuum are evaluated with ONETEP. The solvation energies were also obtained with ONETEP using a minimal parameter implicit solvent model within the self-consistent calculation. Large-scale DFT calculations were also applied within a more theoretically rigorous free energy approach which can, in principle, obtain the full entropic contributions to free energy change. The method performs a mutation from a molecular mechanical (MM) description to an quantum mechanical (QM) description of a system. As a result a QM correction is added to the relative binding free energy obtained from a thermodynamic integration calculation within the MM description. This approach was combined with an electrostatic embedding model within ONETEP and used to calculate the hydration energies of small molecules. As well as the computation of more accurate energies, large-scale DFT calculation compute the electron density of the entire system. Using electron density analysis approaches, such as the Hirshfeld density analysis, in combination with energy decomposition approaches, such as a second order perturbation estimate of natural bond orbital interactions, both qualitative and quantitative understandings can be gained into the contributions of particular chemical functional groups that define protein-ligand interactions. These two approaches where applied to study complexes of the Phosphodiesterase type 5 protein and used to rank ligand binding affinities that agree well with then experimentally observed trends.

Semiconductor and metallic nanoparticles have recently become an attractive area of intensive research due to their unique and diverse properties, that differ significantly from bulk materials. With a wide range of applications and potential uses in nanoelectronics, catalysis, medicine, chemistry or physics an important amount of experimental and theoretical investigations aim to facilitate deeper understating in their physical and chemical behaviour. Within this context, this thesis is focused on the theoretical investigation of silicon, gold and platinum nanoclusters and nanoalloys, in order to provide support for experimental data obtained from collaborating researchers and scientists. Modelled structures of the above nanoparticles were constructed and studied by using a variety of computational tools such as, classical force field MD (DL POLY [1]), tightbinding DFT (DFTB+ [2]), conventional DFT (CASTEP [3]) and linear-scaling DFT (ONETEP [4]). A brief introduction regarding some basic principles of quantum mechanics (QM) and of solid state physics is presented in the first chapter; followed by a general chapter about the classical molecular dynamics (MD) method and its utilisation within the DL POLY code [1]. The last part of the second chapter is devoted to the introduction, validation and implementation of a non-default force field in the source code of DL POLY. The third chapter contains a brief description of some important theorems and terms used in density functional theory (DFT), with some basic information about linear-scaling DFT, as developed in the ONETEP code [4], and tight-binding DFT, reported in the last sections. Chapter 4, includes the results of a series of DFT calculations performed on silicon nanorods, with diameters varying from 0.8 nm to 1.3 nm and about 5.0 nm long. While up to now, similar computational works were conducted on periodic nanowires, in our case, the calculations were performed on the entire nanorods without imposing any symmetry. The fifth chapter proposes a new methodology for calculating extended x-ray absorption fine structure (EXAFS) spectra from modelled geometries of gold nanoparticles by exploiting some of the capabilities of the FEFF code [5]. From several snap-shots of a classical MD simulation, a probability distribution function is calculated for sampling the photoabsorbing and the scattering atoms of the simulated system. The results are then compared with experimental EXAFS data showing a good agreement between the predicted and the measured structures. Finally, in the last two chapters, classical MD simulations on gold and platinum nanoparticles and nanoalloys are reported, which have been performed to support the structural characterisation and analysis of synthesised gold and platinum nanoparticles. Within this framework, DFT calculations have also been attempted on ultrasmall gold nanoparticles and on gold nanosurfaces with one or two thiols attached to them, as a preliminary stage towards the application of linear-scaling DFT in simulating the properties of large metallic systems, currently being studied with semi-empirical quantum approaches or empirical force fields

Kohn-Sham Density Functional Theory (DFT) provides a method for electronic structure calculations applicable to a wide variety of systems. Traditional implementations of DFT are cubic-scaling which limits the size of the systems that can be studied. However recently developed linear-scaling methods, such as onetep, are available which allow much larger systems to be considered. Regardless of scaling DFT has limitations as the exact exchange-correlation functional (a key term in the Kohn-Sham equations) is not known and so approximations have to be made. These approximate functionals generally describe dispersion interactions poorly. In this thesis empirical corrections for dispersion have been developed with parameters optimised for a large set of dispersion bound complexes for the onetep code. This provides a much improved description of dispersion forces which are especially important for biological systems. There is a hierarchy of exchange-correlation functionals available the most accurate of which include a portion of Hartree-Fock exchange. Methods for calculating Hartree- Fock exchange energy in onetep have been developed and are described in this thesis. A quadratic-scaling method using Fourier transforms has been implemented as a benchmark for other implementations. Hartree-Fock exchange may be calculated in a linear-scaling manner by using a numerical pointwise or auxiliary basis set method. Spherical waves have been used as an auxiliary basis set. Linear-scaling has been demonstrated for a polythene chain for these methods. Several hybrid functionals have also been implemented in onetep. These have been validated by comparison with a Gaussian basis set approach in calculations on the reaction paths of an organometallic system