Development and analysis of Tinker-OpenMM as a GPU-based free energy perturbation engine

The utilization of computational technologies for the lead optimization process is one of the biggest challenges in the computational chemistry field. In this dissertation, I describe the addition of GPU-based absolute and relative free energy calculation methods using polarizable force field AMOEBA to Tinker-OpenMM. I then proceed to test the capabilities of this platform by studying the binding free energy and binding structures of derivatives of the MELK inhibitor IN17. Also, I present the implementation of virial-based pressure control to the Tinker-OpenMM platform that is needed for performing isobaric simulations.Cellular and Molecular Biolog

N-body simulations of gravitational dynamics

We describe the astrophysical and numerical basis of N-body simulations, both of collisional stellar systems (dense star clusters and galactic centres) and collisionless stellar dynamics (galaxies and large-scale structure). We explain and discuss the state-of-the-art algorithms used for these quite different regimes, attempt to give a fair critique, and point out possible directions of future improvement and development. We briefly touch upon the history of N-body simulations and their most important results.Comment: invited review (28 pages), to appear in European Physics Journal Plu

Thermodynamic and Structural Phase Behavior of Colloidal and Nanoparticle Systems.

We design and implement a scalable hard particle Monte Carlo simulation toolkit (HPMC), and release it open source. Common thermodynamic ensembles can be run in two dimensional or three dimensional triclinic boxes. We developed an efficient scheme for hard particle pressure measurement based on volume perturbation. We demonstrate the effectiveness of low order virial coefficients in describing the compressibility factor of fluids of hard polyhedra. The second virial coefficient is obtained analytically from particle asphericity and can be used to define an effective sphere with similar low-density behavior. Higher-order virial coefficients --- efficiently calculated with Mayer Sampling Monte Carlo --- are used to define an exponential approximant that exhibits the best known semi-analytic characterization of hard polyhedron fluid state functions. We present a general method for the exact calculation of convex polyhedron diffraction form factors that is more easily applied to common shape data structures than the techniques typically presented in literature. A proof of concept user interface illustrates how a researcher might investigate the role of particle form factor in the diffraction patterns of different particles in known structures. We present a square-triangle dodecagonal quasicrystal (DQC) in a binary mixture of nanocrystals (NCs). We demonstrate how the decoration of the square and triangle tiles naturally gives rise to partial matching rules via symmetry breaking in layers perpendicular to the dodecagonal axis. We analyze the geometry of the experimental tiling and, following the ``cut and project'' theory, lift the square and triangle tiling pattern to four dimensional space to perform phason analysis historically applied only in simulation and atomic systems. Hard particle models are unsuccessful at explaining the stability of the binary nanoparticle super lattice. However, with a simple isotropic soft particle model, we are able to demonstrate seeded growth of the experimental structure in simulation. These simulations indicate that the most important stabilizing properties of the short range structure are the size ratio of the particles and an A--B particle attraction.PhDMaterials Science and EngineeringUniversity of Michigan, Horace H. Rackham School of Graduate Studieshttps://deepblue.lib.umich.edu/bitstream/2027.42/120906/1/eirrgang_1.pd

ACCELERATED COMPUTING FOR MOLECULAR DYNAMICS SIMULATION

Molecular dynamics (MD) simulation serves as a computational microscope into the behavior of the biological and chemical macromolecules. At its core, MD models the interactions between atoms at various levels – force fields model the higher quantum level interactions using simpler physics-based models of interaction energies, while periodic boundary conditions model the bulk phase using lattice-based periodic copies of the simulation box. One limitation of the finite size of the simulation box seen during the simulation of membrane bilayers is the artifact of a chemical disequilibrium between the two layers as a drug molecule enters into the bilayer. We have tried to solve this problem by using a periodic boundary condition which has a half screw symmetry. Our results show that the method scales similar to the best-known method for the normal periodic boundary conditions. We have migrated CHARMM to an efficient implementation on the GPUs. These architectures provide thousands of cores on the same chip but require different programming model in order to use the underlying architecture. Our results show that the new CHARMM CUDA engine is efficient in time and accurate in precision. We have also participated in blind prediction challenges organized by SAMPL community to have a fair assessment of the computational chemistry tools. We developed a hybrid QM and MM technique to predict the pKa of drug-like molecules. It avoids the implicit solvent model used by quantum mechanical models and uses explicit solvent molecules. Since modeling explicit solvent molecules is difficult at QM level, they are modeled at the MM level instead. Thermodynamic cycle couples the aqueous Gibbs free energy of deprotonation to simpler components which can be modeled with higher accuracy. We also built a deep learning model to predict the logP of a set of drug-like molecules in a blind fashion. The generated model is robust over a large number of molecules, not just the ones that it was tested for in the SAMPL competition. We expect the method to be interesting for the drug design industry since lipophilicity of a molecule is important to be known even before it has been synthesized

Quantum mechanical free energy calculations using path integral molecular dynamics

Free energy calculations are one of the most powerful tools within modern theoretical chemistry and are often used to make comparisons with experimental results. Existing free energy calculations are typically performed for classical molecular dynamics simulations but there are certain systems where nuclear quantum effects play an integral role. Specifically, systems with light atoms or low temperatures are the most influenced by such nuclear quantum effects and the development of Feynman path integrals [1] has been effective in accurately describing the quantum nature of these nuclei [2–8]. The primary objective of this thesis is the development of a pair of methodologies to calculate free energies utilizing path integral molecular dynamics to account for nuclear quantum effects. Prior to the development of these free energy methodologies, this thesis presents a communication interface between the OpenMM and MMTK software packages that has been previously published [9]. This interface allows for users of MMTK to take advantage of the performance of OpenMM without major modifications to existing simulation scripts. Notably, the serial OpenMM integrator is shown to provide a 3x performance gain in comparison to a standard MMTK simulation while the GPU implementations of OpenMM provide over a 400x performance gain for larger systems with periodic boundary conditions. The first path integral free energy methodology of this thesis combines the existing um- brella sampling technique [10,11] with path integral molecular dynamics. This methodology has been previously published and proposes that the umbrella sampling biasing potential only needs to be applied to a single path integral bead [12]. Furthermore, this proposed methodology is successfully benchmarked for a pair of Lennard-Jones dimer systems before being applied to the more difficult water dimer. The free energy profiles obtained from simulation are then used to calculate a free energy difference of -12.90 ± 0.05 kJ/mol for the MB-Pol potential in comparison to the experimental dissociation energy of -13.2 ± 0.12 kJ/mol [13]. The second path integral free energy methodology introduces a constraint within the path integral molecular dynamics simulations as opposed to an umbrella sampling restraint. Specifically, this methodology applies a constraint to an individual path integral bead in a manner that is similar to the concept of thermodynamic integration for classical simulations [14]. Formal estimators for the derivative of the free energy have been developed by Iouchtchenko et al. [15] and the results presented in this thesis analyze the effectiveness of these estimators for molecular dynamics simulations of Lennard-Jones and water dimers. Additionally, a new estimator is developed and the resulting free energy profiles are used to evaluate a free energy difference for the water dimer of -13.03 ± 0.14 kJ/mol, which is within the errors of the experimental dissociation energy [13]. Overall, this thesis provides a theoretical framework to study the free energy of weakly bound systems over a broad range of temperatures. It is important to note that these methodologies were insufficient below 25 K and it remains more practical to use reaction coordinates that are not distances at such temperatures. Nevertheless, the extension and application of these methodologies to more complicated systems remains an area of exciting development

Development of High Performance Molecular Dynamics with Application to Multimillion-Atom Biomass Simulations

An understanding of the recalcitrance of plant biomass is important for efficient economic production of biofuel. Lignins are hydrophobic, branched polymers and form a residual barrier to effective hydrolysis of lignocellulosic biomass. Understanding lignin\u27s structure, dynamics and its interaction and binding to cellulose will help with finding more efficient ways to reduce its contribution to the recalcitrance. Molecular dynamics (MD) using the GROMACS software is employed to study these properties in atomic detail. Studying complex, realistic models of pretreated plant cell walls, requires simulations significantly larger than was possible before. The most challenging part of such large simulations is the computation of the electrostatic interaction. As a solution, the reaction-field (RF) method has been shown to give accurate results for lignocellulose systems, as well as good computational efficiency on leadership class supercomputers. The particle-mesh Ewald method has been improved by implementing 2D decomposition and thread level parallelization for molecules not accurately modeled by RF. Other scaling limiting computational components, such as the load balancing and memory requirements, were identified and addressed to allow such large scale simulations for the first time. This work was done with the help of modern software engineering principles, including code-review, continuous integration, and integrated development environments. These methods were adapted to the special requirements for scientific codes. Multiple simulations of lignocellulose were performed. The simulation presented primarily, explains the temperature-dependent structure and dynamics of individual softwood lignin polymers in aqueous solution. With decreasing temperature, the lignins are found to transition from mobile, extended to glassy, compact states. The low-temperature collapse is thermodynamically driven by the increase of the translational entropy and density fluctuations of water molecules removed from the hydration shell

CP2K: An electronic structure and molecular dynamics software package - Quickstep: Efficient and accurate electronic structure calculations

CP2K is an open source electronic structure and molecular dynamics software package to perform atomistic simulations of solid-state, liquid, molecular, and biological systems. It is especially aimed at massively parallel and linear-scaling electronic structure methods and state-of-the-art ab initio molecular dynamics simulations. Excellent performance for electronic structure calculations is achieved using novel algorithms implemented for modern high-performance computing systems. This review revisits the main capabilities of CP2K to perform efficient and accurate electronic structure simulations. The emphasis is put on density functional theory and multiple post–Hartree–Fock methods using the Gaussian and plane wave approach and its augmented all-electron extension

Role of the anisotropy in the interactions between nano- and micro-sized particles

The present Thesis focuses on the thermodynamic and dynamic behaviour of anisotropically interacting colloids by means of theoretical and numerical techniques. Colloidal suspensions, i.e. micro-- and nano--sized particles dispersed in a continuous phase, are a topic of great interest in several fields, including material science, soft matter and biophysics. Common in everyday life in the form of soap, milk, cream, etc., colloids have been used for decades as models for atomic and molecular systems, since both classes of systems share many features like critical phenomena, crystallisation and glass transition. Experimental investigation of colloidal systems is made easier by the large size of colloids, which makes it possible to employ visible light as an experimental probe to investigate these systems. Moreover, since the mass of the particles controls the timescales of the dynamics, relaxation times of colloidal suspensions, ranging from seconds to years, orders of magnitude larger than their atomic counterparts, are more easily experimentally accessible. By exploiting this intrinsic slowness, with respect to molecular liquids, present day experimental techniques make it possible to follow in time trajectories of ensembles of particles with tools like confocal microscopy, thus effectively allowing to reconstruct the whole phase space trajectory of the system. In addition, it is also possible to manipulate single and multiple objects using techniques like optical tweezers, magnetic tweezers and atomic force microscopy. With single-molecule force spectroscopy one can arrange particles in ordered structures or measure properties like stiffness or mechanical responses (as in pulling experiments on RNA and DNA strands of particles and aggregates). A remarkable difference between the molecular and the colloidal world is that in the former the interactions between the basic constituents are fixed by nature, while in the latter the effective potential between two particles can be controlled by accurately designing and synthesizing the building blocks or tuned by changing the properties of the solvent. In the last decade many new sophisticated techniques for particle synthesis have been developed and refined. These recent advances allow for the creation of an incredible variety of non-spherically, i.e. anisotropically, interacting building blocks. The anisotropy can arise from shape, surface patterning, form of the interactions or a combination thereof. Examples are colloidal cubes, Janus particles, triblock Janus particles, patchy particles, magnetic spheres and many others. The recent blossoming of experimental, theoretical and numerical studies and research on the role of the anisotropy has highlighted the richness of phenomena that these systems exhibit. Relevant examples for the present Thesis are valence-limited building blocks, i.e colloids with a maximum number of bound neighbours, and non-spherical particles with an aspect ratio, i.e. the ratio of the width of a particle to its height, significantly different from $1$. The simplest example of valence-limited colloids is given by the so-called \textit{patchy} particles: colloids decorated with attractive spots (patches) on the surface. If the width and the range of the patches are chosen in such a way that each patch can form no more than one bond, then the total number of bound first neighbours per particle $M$ can not exceed the number of patches. For particles interacting through short-ranged isotropic potentials, $M \approx 12$. It has been shown that changing the valence $M$ has dramatic effects, both qualitative and quantitative, on the dynamic and thermodynamic properties of such systems. At high densities patchy colloids can self-assemble into a large variety of crystal structures, depending on valence, geometry and external parameters. We will mostly focus on low-density systems. The second class of systems pertinent to the present work comprises anisotropically shaped particles that, depending on the aspect ratio and the values of the external parameters, can exhibit liquid crystal phases which may display orientational long-range order. Nematic, in which there is no translational order, smectic, in which particles are ordered in layers and thus exhibit translational order in one dimension, and columnar phases, in which particles self-assemble into cylindrical aggregates which can in turn become nematic or form two-dimensional lattices, do not exist in isotropic systems, since the anisotropy in shape is a prerequisite for the breaking of the orientational symmetry. Liquid crystals, discovered at the end of the 19th century have been thoroughly investigated for decades, leading to technological breakthroughs like LCD displays. Recently it has been suggested that liquid crystal phases occurring in dense solutions of short DNA double strands could have played a role in the prebiotic chemical generation of complementary H-bonded molecular assemblies. The main goal of the present Thesis is to study the structural, thermodynamic and, to a lesser extent, dynamic properties of systems interacting through anisotropic potentials at low densities and temperatures. In particular, we focus on the low-density phase behaviour of valence-limited systems. We use a variegated approach, comprising state-of-the-art Monte Carlo and Molecular Dynamics techniques and theoretical approaches, to analyse and shed some light on the effect of the anisotropy on the phase diagram and on the dynamics of such systems. As the effect of the valence on the phase diagram plays a major role in the models investigated throughout this Thesis, each Chapter is devoted to the study of the dynamics and thermodynamics of systems having a fixed or effective maximum valence $M$. In the last years a lot of effort has been devoted to the study of end-to-end stacking interactions between different strands of nucleic acids, which play an important role in both physical and biological applications of DNA and RNA. In Chapter~1, building on the experimental work of Bellini \textit{et al.}, we make use of a theoretical framework recently developed to tackle the problem of the isotropic--nematic phase coexistence in solutions of short DNA duplexes (DNADs). We compare the parameter-free theoretical predictions with results from large scale numerical simulations on GPUs of a coarse-grained realistic model and find a good quantitative agreement at low concentrations. We then predict the phase boundaries for different DNAD lengths and compare the results with experimental findings. In Chapter~2 we investigate the structural and thermodynamic properties of systems having $M=2$, that is systems that undergo an extensive formation of linear structures as temperature is lowered. We focus on bi-functional patchy particles whose interaction details are chosen to qualitatively mimic the behaviour of the low-density, low-temperature dipolar hard sphere (DHS) model by analysing the outcomes of the simulations carried out in Chapter~3. In particular, we are interested in the interplay between chains and rings in equilibrium polymerization processes in a region of the phase diagram where the formation of the latter is favoured. The very good quantitative agreement found by comparing numerical results with theoretical, parameter-free predictions calls for an extension of the theory with the inclusion of branching, in order to understand how the presence of rings affects the phase separation. Chapter~3 is devoted to the investigation of the phase behaviour of dipolar fluids, i.e. systems interacting mainly through dipole-dipole potentials. For spheres, the lowest-energy configuration is the nose-to-tail contact geometry, and hence the ground state is an infinite chain or ring like in regular $M=2$ systems. For finite temperatures, on the other hand, thermal fluctuations allow for the appearance of defects like dangling ends and chain branching which, in the language of this Thesis, makes for a temperature-dependent valence. This general mechanism, under some specific conditions, can lead to a very peculiar phase separation, driven by a balance between these \textit{topological} defects rather than by the energy/entropy competition usually responsible for regular gas--liquid phase transitions. This topological phase transition has been recently observed in a model system of patchy particles but it is unclear whether such mechanism still holds in dipolar fluids in general and in the DHS model in particular. We focus on the DHS model, whose phase behaviour at low densities and temperatures has been studied for decades but still remains largely unknown. In particular, we look for the gas--liquid critical point by means of state-of-the-art Monte Carlo simulations in a region where it has long been thought to be. We find no evidence of a phase transition and we speculate that this is due to an abundance of rings, providing a remarkable example of phase separation suppressed by self-assembly. In Chapter~4 we study the dynamics of tetravalent patchy particles in the optimal network density region. For this fixed value of density the system is able to form a fully connected random network, i.e. an ideal gel. Indeed, as the temperature is lowered, a percolating network forms and the dynamics slows down. Although the observed dynamical arrest is different from the glass case, where excluded volume interactions are dominant, the decay of the self-- and collective correlation functions of the resulting fluid bears similarities with that observed in glassy systems. Remarkably, comparing the characteristic decay times of density-density correlation functions with the average bond life, we find that only at very low $T$ the decay of the density fluctuations requires the breakage of bonds. In Chapter~5 we introduce DNA as a building block that can be used to rationally design novel, self-assembling materials with tunable properties. In this Chapter, we study the phase behaviour and the dynamics of four-armed DNA constructs at low densities. We use the coarse-grained, realistic DNA model employed in Chapter~1 and state-of-the-art simulation techniques, as presented in Chapter~6, to investigate systems composed of thousands of nucleotides undergoing a two-step self-assembling process and we quantitatively compare the outcome with experimental results obtained for a very similar system. In Chapter~6 we introduce Graphics Processing Units (GPUs) as valuable tools for present day numerical investigations. We outline both the architecture of NVIDIA GPUs and NVIDIA CUDA, the software layer built on top of the hardware required to program these devices. We then present the techniques employed to write an efficient, general Molecular Dynamics code and compare its performances with a regular CPU code. The observed performance boost allows us to tackle the analysis of the dynamics and thermodynamics of very large systems without having to resort to massive CPU clusters (see Chapters~1,~4 and~5). Our work shows that it is possible to predict the location of thermodynamic and dynamic \textit{locii} of very complicated objects by means of numerical simulations. Since the available computational power keeps increasing at a steady pace, it will be soon possible to repeat the pioneering study presented in this Thesis on a more automated basis and for even more complicated system. For example, it will be possible to directly study the isotropic--nematic phase transition of short DNA duplexes investigated in Chapter~1 or design self-assembling DNA strands able to reproduce the behaviour of the patchy colloids or dipolar fluids studied throughout this Thesis. Being able to carefully design the building blocks and then predict beforehand the properties of a compound will greatly simplify the process of synthesising tomorrow's materials

Nuclear Quantum Effects in Water and Aqueous Systems: Experiment, Theory, and Current Challenges

Nuclear quantum effects influence the structure and dynamics of hydrogen-bonded systems, such as water, which impacts their observed properties with widely varying magnitudes. This review highlights the recent significant developments in the experiment, theory, and simulation of nuclear quantum effects in water. Novel experimental techniques, such as deep inelastic neutron scattering, now provide a detailed view of the role of nuclear quantum effects in water's properties. These have been combined with theoretical developments such as the introduction of the principle of competing quantum effects that allows the subtle interplay of water's quantum effects and their manifestation in experimental observables to be explained. We discuss how this principle has recently been used to explain the apparent dichotomy in water's isotope effects, which can range from very large to almost nonexistent depending on the property and conditions. We then review the latest major developments in simulation algorithms and theory that have enabled the efficient inclusion of nuclear quantum effects in molecular simulations, permitting their combination with on-the-fly evaluation of the potential energy surface using electronic structure theory. Finally, we identify current challenges and future opportunities in this area of research